The problem can be solved using only elementary arguments (i.e without RMT).
Claim. In the limit $d,m \to \infty$ such that $m/d \to \rho \in (0,\infty)$, it holds that $$ m^{-1}d^{-2}\sum_{i}\lambda_i(S)^2 \overset{a.s}{\to} 1+\rho. $$
Indeed, one may write $$ \sum_{i}\lambda_i(S)^2 = \sum_{i,j=1}^m s_{i,j}^2 = \sum_{i,j} (a_i^\top a_j)^2 = \sum_{i=1}^m(\|a_i\|^4 + \sum_{j\ne i}^m (a_i^\top b_j)^2). \tag{1} $$ Now, by the Law of Large numbers, it's clear that $d^{-2}\|a_i\|^4 = (d^{-1}\|a_i\|^2)^2 \overset{a.s}{\to} 1^2 = 1$ for all $i \in [m]$. On the other hand, if $i \ne j$, then $d^{-2}(a_i^\top b_j)^2$ has a beta distribution with parameters $\alpha=1$ and $\beta=d-2$ (see this post https://mathoverflow.net/a/227156/78539), and expected value $\alpha / (\alpha + \beta) = 1/(d-1)$. Combining with (1) via another application of LLNs (to handle the second term) then completes the proof after noting that $(m-1) / (d-1) \to \rho$.
Question. Comparing with the accepted answer, is true that $\langle \lambda^2\rangle_{MP(1/\rho)} = 1 + \rho$ ?