What is the relation between the eigenvectors of a sample covariance matrix and those of the true covariance matrix?

As is known, the covariance matrix of a set of random vectors $\{\mathbf{x}_i\}_{i=1}^N$ can be estimated by their sample covariance matrix: $\mathbf{\hat R}:=\frac{1}{N}\sum_{n=1}^N\mathbf{x}_n\mathbf{x}_n^T$
In popular subspace methods such as PCA and MUSIC, the eigenvectors of $\mathbf{\hat R}$ is used to estimate the signal subspace and noise subspace.
As far as I know, there are mainly two streams of results about the relation between the sampled eigenvectors $\mathbf{\hat U}:=[\mathbf{\hat u}_1,\mathbf{\hat u}_2,\cdots]$ of $\mathbf{\hat R}$ and the true eigenvectors $\mathbf{U}:=[\mathbf{u}_1,\mathbf{u}_2,\cdots]$ of $\mathbf{R}$,
Firstly, under Gaussian assumption, the sampled eigenvectors $\mathbf{\hat U}$ obey an asymptotic normal distribution (in the sense of large samples N)
$$\mathbf{\hat U}\simeq \mathcal{N}(\mathbf{U},\mathbf{\Sigma})$$
these results were conducted by T.W. Anderson, P. Stoica, etc.
And secondly, the sampled eigenvectors $\mathbf{\hat U}$ can be modeled as a perturbation version of the true eigenvectors $\mathbf{U}$, that is $$\mathbf{\hat U}=\mathbf{U} + \delta\mathbf{U}$$
where the perturbation $\delta\mathbf{U}$ is related with the perturbation caused by finite samples, these results were conducted by G Stewart, F Li, etc.
My question is, how should we evaluate these two kinds of analysis, I mean, which one is better from what aspect, and why?
And, are there new and more advanced results about the relationship between the $\mathbf{\hat U}$ and $\mathbf{U}$ mentioned above? Is there any general theoretical background lying behind these analyzes?

The first focuses on what happens when the number of samples tends to infinity, and focuses on the distribution of the eigenvectors (assuming a specific underlying distribution of the $\boldsymbol{x}_i$). It is a probabilistic bound; you could get any estimated covariance matrix $\tilde{R}$ from the observations, and any set of eigenvectors, but some are more likely than others.
The second bound focuses on how the eigenvalues react to a perturbation in $R$. It tells you how you should expect the eigenvalues to change if you have a different matrix $\tilde{R}$ rather than $R$. This different matrix might come from the estimation process, but the bound has nothing probabilistic in nature. We have one matrix $\tilde{R}$, close to $R$, and we are interested in how distant and in which direction its eigenvectors are with respect to those of $R$.
If anything, the second bound can be useful in proving the first one, as the final part of the proof: first you prove that the sample covariance $\tilde{R}$ is asymptotically normal, then apply a delta method-type argument, using the fact that the eigenvectors are a linear (in first-order) function of the perturbation (which is essentially what the Stewart-Li-type bounds state, under the assumption that the eigenvalues are simple).