For $N,n\in\mathbb N$, consider the sample covariance matrix $$S_N(T)=\frac{1}{N}XX^*$$ from a data set $X=C_N^{1/2}ZT_n^{1/2}$, where $Z=(Z_{i,j})$ is a $N\times n$ matrix having i.i.d. entries with mean zero and variance one, and $C_N, T_n$ are deterministic positive semi-definite Hermitian matrices of dimension $N$ and $n$, respectively. We assume that $(C_N)_N$ is bounded in spectral norm, and $T_n$ is a Toeplitz matrix with its largest eigenvalues diverging to infinity. The matrix $X$ can be viewed as a data set of an $N$-dimensional long memory stationary process having separable dependence structure. As $N,n\to \infty$ and $Nn^{-1} \to r\in (0,\infty)$, we establish the asymptotics and the joint CLT for $(λ_1(S_N(T)),\cdots, λ_m(S_N(T))^\top$ where $λ_j(S_N(T))$ denotes the $j$th largest eigenvalue of $S_N(T)$, and $m$ is a fixed integer. For the CLT, we first study the case where the entries of $Z$ are Gaussian, and then we generalize the result to some more generic cases. This result substantially extends our previous result in Merlevède et al. 2019, where we studied $λ_1(S_N(T))$ in the case where $m=1$ and $X=ZT_n^{1/2}$ with $Z$ having Gaussian entries. In order to establish this CLT, we are led to study the first order asymptotics of the largest eigenvalues and the associated eigenvectors of some deterministic Toeplitz matrices. We are specially interested in the autocovariance matrices of long memory stationary processes. We prove multiple spectral gap properties for the largest eigenvalues and a delocalization property for their associated eigenvectors.