Consider the vector space $\mathbb{R}^n$, the standard inner product $\langle \cdot,\cdot \rangle:\mathbb{R}^n\times\mathbb{R}^n\rightarrow \mathbb{R}$, and some $0<\epsilon\leq \frac{1}{\sqrt{n}}$. Is it possible to generate a set of $M$ vectors, say $\mathcal{S}$, in a deterministic fashion (not random), such that the vectors in $\mathcal{S}$ satisfy the following properties:

all the entries of the vectors come from $\{\frac{-1}{\sqrt{n}},\frac{1}{\sqrt{n}}\}$, so that $||\mathbb{x}||_2=1\,\,\forall \mathbb{x}\in\mathcal{S}$,

$\mathbb{x},\mathbb{y}\in\mathcal{S}\Rightarrow|\langle \mathbb{x},\mathbb{y} \rangle| <\epsilon$, i.e., the vectors in $\mathcal{S}$ are almost orthogonal to each other.

Kindly note that a similar question was asked before (Almost orthogonal vectors in subsets of euclidean space), but it asked for an upper bound on $M$ for given $n,\epsilon$ - the answer was $M\leq \frac{n(1-\epsilon^2)}{1-n\epsilon^2}$. I am interested in knowing if there is a deterministic way of constructing $\mathcal{S}$ of reasonable size.

Let me also explain an approach that I thought of (but got stuck) - consider $\mathcal{C}=[n, k, d]_{2}$, a binary linear code with $\frac{n}{2}(1-\epsilon)\leq d <\frac{n}{2}(1+\epsilon)$, and map its codewords (of length $n$) to real-valued sequences in $\{\frac{-1}{\sqrt{n}},\frac{1}{\sqrt{n}}\}^n$ such that $0$ maps to $\frac{-1}{\sqrt{n}}$ and $1$ maps to $\frac{1}{\sqrt{n}}$. Then consider the set of real-valued sequences, say $\mathcal{T}$, corresponding to the subset of codewords in $\mathcal{C}$ whose Hamming weight is at least $d$ and at most $\frac{n}{2}(1+\epsilon)$. Then one can show that $\mathcal{T}$ is a candidate for $\mathcal{S}$; however I have not been able to design a linear code for which tight upper and lower bounds on $|\mathcal{T}|$ can be obtained (i.e., I don't have bounds/expression for $M$ in terms of $n,\epsilon$).