Let $X,X_1,\dots, X_n$ be i.i.d. Gaussian random variables in a separable Hilbert space ${\mathbb H}$ with zero mean and covariance operator $Σ={\mathbb E}(X\otimes X),$ and let $\hat Σ:=n^{-1}\sum_{j=1}^n (X_j\otimes X_j)$ be the sample (empirical) covariance operator based on $(X_1,\dots, X_n).$ Denote by $P_r$ the spectral projector of $Σ$ corresponding to its $r$-th eigenvalue $μ_r$ and by $\hat P_r$ the empirical counterpart of $P_r.$ The main goal of the paper is to obtain tight bounds on $$ \sup_{x\in {\mathbb R}} \left|{\mathbb P}\left\{\frac{\|\hat P_r-P_r\|_2^2-{\mathbb E}\|\hat P_r-P_r\|_2^2}{{\rm Var}^{1/2}(\|\hat P_r-P_r\|_2^2)}\leq x\right\}-Φ(x)\right|, $$ where $\|\cdot\|_2$ denotes the Hilbert--Schmidt norm and $Φ$ is the standard normal distribution function. Such accuracy of normal approximation of the distribution of squared Hilbert--Schmidt error is characterized in terms of so called effective rank of $Σ$ defined as ${\bf r}(Σ)=\frac{{\rm tr}(Σ)}{\|Σ\|_{\infty}},$ where ${\rm tr}(Σ)$ is the trace of $Σ$ and $\|Σ\|_{\infty}$ is its operator norm, as well as another parameter characterizing the size of ${\rm Var}(\|\hat P_r-P_r\|_2^2).$ Other results include non-asymptotic bounds and asymptotic representations for the mean squared Hilbert--Schmidt norm error ${\mathbb E}\|\hat P_r-P_r\|_2^2$ and the variance ${\rm Var}(\|\hat P_r-P_r\|_2^2),$ and concentration inequalities for $\|\hat P_r-P_r\|_2^2$ around its expectation.