A possible relevant post What kind of random matrices have rapidly decaying singular values?. In that post I discussed the distribution of maximal eigenvalue of a random matrix based on the result [Johnstone]. However that answer does not fully answer your question since there I emphasized the fact that the maximal eigenvalues of most random matrices in applications(Wishart and Uniform on a Grassmanian) will follow Tracy-Widom distribution asymptotically as long as the size of matrix grows with the $p$ in a reasonable manner.
$\blacksquare$1. Why the concentration bound does not work for Rademacher example?
This post posed a slightly more complicated problem. $\bar X_{n}$ does not necessarily be Hermitian nor almost bounded in each of its entries. As you observed that assumptions are always put onto the random matrices in order to make it into a covariance matrix from some random vector.([Mackey] and [Tropp] requires Hermitian in order to use Stein-pair argument;[Vershynin] mostly requires a lot of isotropic conditions to reduce his argument to a lower dimension where the tail is tractable.) Most of the existing results about the concentration bound focus on the case where the random matrix comes from a covariance matrix because of not only applied interest but also its positivity. See a discussion in [Guédon et al.].
This is the first obstacle why the results you pointed out might not be applicable to $\lambda_{max}(\bar{X})$. The example you constructed using Rademacher random variables has the problem that it is not positive definite, in fact its rank is always one regardless of the dimension of $x_i$ you choosen, therefore it will correspond to a degenerate distribution whose support will concentrate on a lower dimension space.
I do not think your comparison is of any meaning by "applying" results for covariance matrices onto your example or Rademacher r.v.s. In fact in your example, suppose we have $x_i=(r_{i1},r_{i2},r_{i3})$ where $r_{ij}$ are iid Rademacher random variables, then it is readily seen that the eigenvalues of $X_i=x_i x_i^{T}$ will always be $(3,0,0)$, in that case it is meaningless to talk about converging rate. For $x_i=(r_{i1},r_{i2},\cdots, r_{in})$ it is not hard to verify it has eigenvalues of $x_i x_i^{T}$ will be $n$ with multiplicity one and $0$ with multiplicity $(n-1)$. So no matter how larger the Rademacher r.v. is, it is actually a one dimensional degenerated distribution and $X_i=x_i x_i^{T}$ cannot be a covariance matrix.
Therefore to your question
My question is why the bounds obtained in random matrix theory
literature do not yield the sharpest convergence rates for these
simple cases? Is this a by-product of the employed techniques or is
there something fundamental that prevents it ?
The semi-positivity of "simple cases" that prevents it from getting the sharpest rate control. As for the case where $x_i$ are uniformly distributed random variables with zero mean, we can also prove that its eigenvalues are $O(1)$(big O) with multiplicity 1 and $o(1)$(small O) with multiplicity $(n-1)$ so that positivity is preserved yet only one eigenvalue contribute to the concentration bound in $O(1)$ order.
So I think the failure of attaining sharpest bound is due to the nature of your example instead of the technique applied in the proof.
$\blacksquare$2. What is the concentration bound for Gaussian matrices?
A sharp concentration bound is available in [Bandeira&van Handel].
Back to the second part of your question, you are correct that $\lambda_{max}(X_i)≤R$ will not even allow Gaussian random matrices, but we usually study a more general statistical model called matrix model[Pastur et.al], which has densities in form of $$P_{n,\beta}(d_{\beta}M)=Z_{n,\beta}^{-1}exp\left\{ -\frac{1}{2}n\beta trV(M)\right\} d_{\beta}M$$[Pastur et.al] (4.1.1) where $V(\lambda)\geq(2+\epsilon)log(1+|\lambda|)$ is a continuous function ensuring the integrability of the matrix measure and where $\beta=1$ is equivalent to the real symmetric matrices and $\beta=2$ is equivalent to Hermitians. For the Gaussian model you pointed ou, we can realized it for $V(\lambda)=\frac{\lambda^2}{2w}$ so if you admit that integrability is essential to such a matrix model then the concentration bound must be of the form shown and cannot even be discrete.[Pastur et.al] pp.105-106
In fact the eigenvalues for Gaussian matrices has following distribution density as derived in Pastur et.al$$\tilde{Q}_{n,\beta}^{-1}exp\left\{ -\frac{n\beta}{2}\sum_{l=1}^{n}V(\lambda_{l})\right\} \left|\prod_{1\leq j<k\leq n}\left(\lambda_{j}-\lambda_{k}\right)\right|^{\beta}\prod_{j=1}^{n}d\lambda_{j}\left[H_{\beta}\right]$$
So the bound must not even be sharper than this distribution law, which is providing an exponential bound.
If as you said in your post, the $X_i$'s are self-adjoint, then in Gaussian case the $\bar{X}$ will again be distributed as a Gaussian distribution, and the distribution above will also apply to your $\bar{X}$. Generally, [Johnstone] gave a relatively good bound (under a mild constraint on dimension $p$). See the linked post (What kind of random matrices have rapidly decaying singular values?).
My point here is that we should first figure out the distribution of $\bar{X}$ and then discuss what kind of concentration bound applies. It is also possible that although the assumptions of concentration bound are satisfied by $X_i$ but not satisfied by $\bar{X}$
$\blacksquare$3.Extension of [Kutin]'s weakly bounded difference result.
After a careful examination of [Kutin], I found your argument not very convincing. The weakly difference bounded conditions (Def. 2.3 in [Kutin]) also required that $|X(\omega)-X(\omega')|<b_k$ when $\omega$ and $\omega'$ differ only in the k-th coordinate besides the McDiarmid-type condition $$Pr_{(\omega,v)\in\Omega\times\Omega_{k}}(\left|X(\omega)-X(\omega')\right|<c_{k})\geq1-\delta$$, if you use $\Omega_k$ as the probability spaces associated with the Gaussian random variables $X_k$, although the $\delta$ probability part of the weakly bounded difference condition holds, the exact bounded part $|X(\omega)-X(\omega')|<b_k$ of the weakly bounded condition fails for two independent Gaussian random variable for any bounded sequence ${b_k}$.
So I do not think you can apply [Kutin]'s result directly. And by the way if you have a chance to read this, could you add the literature you mentioned in OP? Thanks
Reference
[Vershynin]http://www-personal.umich.edu/~Vershyninv/papers/non-asymptotic-rmt-plain.pdf
[Tropp]Tropp, Joel A. "User-friendly tail bounds for sums of random matrices." Foundations of computational mathematics 12.4 (2012): 389-434.
[Johnstone]Johnstone, Iain M. "Multivariate analysis and Jacobi ensembles: Largest eigenvalue, Tracy–Widom limits and rates of convergence." Annals of statistics 36.6 (2008): 2638.(Slides: http://www.stat.harvard.edu/NESS10/Johnstone-slides.pdf)
[Pastur et.al]Pastur, Leonid Andreevich, Mariya Shcherbina, and Mariya Shcherbina. Eigenvalue distribution of large random matrices. Vol. 171. Providence, RI: American Mathematical Society, 2011.
[Kutin]Kutin, Samuel. "Extensions to McDiarmid’s inequality when differences are bounded with high probability." Department Computer Science, University of Chicago, Chicago, IL. Technical report TR-2002-04 (2002).
[Bandeira&van Handel]Bandeira, Afonso S., and Ramon van Handel. "Sharp nonasymptotic bounds on the norm of random matrices with independent entries." The Annals of Probability 44.4 (2016): 2479-2506.
[Guédon et al.]Guédon, Olivier, et al. "The central limit theorem for linear eigenvalue statistics of the sum of independent random matrices of rank one." Spectral Theory and Differential Equations. Amer. Math. Soc. Transl. Ser 2.233 (2014): 145-164.