  In this paper, we prove the edge universality of largest eigenvalues for
separable covariance matrices of the form $\mathcal Q :=A^{1/2}XBX^*A^{1/2}$.
Here $X=(x_{ij})$ is an $n\times N$ random matrix with $x_{ij}=N^{-1/2}q_{ij}$,
where $q_{ij}$ are $i.i.d.$ random variables with zero mean and unit variance,
and $A$ and $B$ are respectively $n \times n$ and $N\times N$ deterministic
non-negative definite symmetric (or Hermitian) matrices. We consider the
high-dimensional case, i.e. ${n}/{N}\to d \in (0, \infty)$ as $N\to \infty$.
Assuming $\mathbb E q_{ij}^3=0$ and some mild conditions on $A$ and $B$, we
prove that the limiting distribution of the largest eigenvalue of $\mathcal Q$
coincide with that of the corresponding Gaussian ensemble (i.e. the $\mathcal
Q$ with $X$ being an $i.i.d.$ Gaussian matrix) as long as we have $\lim_{s
\rightarrow \infty}s^4 \mathbb{P}(\vert q_{ij} \vert \geq s)=0$, which is a
sharp moment condition for edge universality. If we take $B=I$, then $\mathcal
Q$ becomes the normal sample covariance matrix and the edge universality holds
true without the vanishing third moment condition. So far, this is the
strongest edge universality result for sample covariance matrices with
correlated data (i.e. non-diagonal $A$) and heavy tails, which improves the
previous results in \cite{BPZ1,LS} (assuming high moments and diagonal $A$),
\cite{Anisotropic} (assuming high moments) and \cite{DY} (assuming diagonal
$A$).
