For N, n is an element of N, consider the sample covariance matrix S-N(T)=1/NXX*. from a data set X = C(N)(1/2)ZT(n)(1/2), where Z = (Z(i,j)) is a N x 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 -> infinity and Nn(-1) -> r is an element of (0, infinity), we establish the asymptotics and the joint CLT for (lambda(1) (S-N (T)), ..., lambda(m)(S-N(T))(inverted perpendicular), where lambda(j)(S-N(T)) denotes the jth largest eigenvalue of S-N(T) and in 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 Merlevede et al. 2019, where we studied lambda(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 related to long memory stationary processes. We prove multiple spectral gap properties for the largest eigenvalues and a delocalization property for their associated eigenvectors.