Let $X,X_1,\dots, X_n$ be i.i.d. Gaussian random variables with zero mean and covariance operator $Σ={\mathbb E}(X\otimes X)$ taking values in a separable Hilbert space ${\mathbb H}.$ Let $$ {\bf r}(Σ):=\frac{{\rm tr}(Σ)}{\|Σ\|_{\infty}} $$ be the effective rank of $Σ,$ ${\rm tr}(Σ)$ being the trace of $Σ$ and $\|Σ\|_{\infty}$ being its operator norm. Let $$\hat Σ_n:=n^{-1}\sum_{j=1}^n (X_j\otimes X_j)$$ be the sample (empirical) covariance operator based on $(X_1,\dots, X_n).$ The paper deals with a problem of estimation of spectral projectors of the covariance operator $Σ$ by their empirical counterparts, the spectral projectors of $\hat Σ_n$ (empirical spectral projectors). The focus is on the problems where both the sample size $n$ and the effective rank ${\bf r}(Σ)$ are large. This framework includes and generalizes well known high-dimensional spiked covariance models. Given a spectral projector $P_r$ corresponding to an eigenvalue $μ_r$ of covariance operator $Σ$ and its empirical counterpart $\hat P_r,$ we derive sharp concentration bounds for bilinear forms of empirical spectral projector $\hat P_r$ in terms of sample size $n$ and effective dimension ${\bf r}(Σ).$ Building upon these concentration bounds, we prove the asymptotic normality of bilinear forms of random operators $\hat P_r -{\mathbb E}\hat P_r$ under the assumptions that $n\to \infty$ and ${\bf r}(Σ)=o(n).$ In a special case of eigenvalues of multiplicity one, these results are rephrased as concentration bounds and asymptotic normality for linear forms of empirical eigenvectors. Other results include bounds on the bias ${\mathbb E}\hat P_r-P_r$ and a method of bias reduction as well as a discussion of possible applications to statistical inference in high-dimensional principal component analysis.