I'm wondering how $\epsilon_{ijk}\epsilon_{ilm} $ can give this matrix $\begin{pmatrix} \delta_{ii} & \delta_{ij} & \delta_{ik} \\ \delta_{li} & \delta_{lj} & \delta_{lk} \\ \delta_{mi} & \delta_{mj} & \delta_{mk} \\ \end{pmatrix}$
Is it a cross product between 2 vectors?
The "Matrix" you are referring to is actually a determinant.
$$\left| \begin{array}{ccc} \delta^i_i & \delta^i_j & \delta^i_k \\ \delta^l_i &\delta^l_j & \delta^l_k \\ \delta^m_i &\delta^m_j & \delta^m_k \\ \end{array} \right|=\left| \begin{array}{ccc} \delta^l_j & \delta^l_k \\ \delta^m_j & \delta^m_k \\ \end{array} \right|=\delta^l_j\delta^m_k-\delta^l_k\delta^m_j$$
If you look closer you will find that the the first index of $\delta$, indices $i,l,m$, correspond to the 1st, 2nd and 3rd rows of this determinant, whereas the second index, $i,j,k$, correspond to the 1st, 2nd and 3rd columns of the determinant.
From introductory linear algebra we remember that an interchange of two rows (or two columns) changes the sign of the determinant; in particular, if two rows (or two columns) are identical, the value of the determinant is zero.
Clearly the determinant is skew-symmetric in both sets of indices (rows $i,l,m$ and columns $i,j,k$). For convenience we can denote this determinant $\delta^{ilm}_{ijk}$.
Now lets study the difference (if any) between the two entities
$$D^{ilm}_{ijk}=\varepsilon^{ilm}\varepsilon_{ijk}-\delta^{ilm}_{ijk}\tag{1}$$
Since the Levi-Civita symbols are also totally skew-symmetric the only nonzero components of $D^{ilm}_{ijk}$ will occur when the index sets $ilm$ and $ijk$ are each a permutation of $1,2,3$.
However, upon inspection of $(1)$ it is easily seen that (when the first term is 1 the second is -1 and vice versa due to the skew symmetric property)
$D^{ilm}_{ijk}\equiv 0$
So we conclude that $$\delta^{ilm}_{ijk}=\varepsilon^{ilm}\varepsilon_{ijk}$$ This line of reasoning is also valid for the more general case
$$\delta^{j_1\dots j_n}_{h_1\dots h_n}=\varepsilon_{h_1\dots h_n}\varepsilon^{j_1\dots j_n}$$
It turns out that the determinant (usually called the generalized Kronecker delta) is powerful tensor in its own right. Unlike the Levi-Civita symbols it can have an arbitrary number of indices and is often extremely useful in complicated manipulative processes.
We should also mention that in the specific case above there is a contraction taking place over index $i$ (due to the summation convention)
$$\varepsilon^{ilm}\varepsilon_{ijk}=\delta^{ilm}_{ijk}=\delta^{lm}_{jk}=\delta^l_j\delta^m_k-\delta^l_k\delta^m_j $$
Where we have used the convenient and beautiful rule for the contraction of a generalized Kronecker delta:
$$\delta^{j_1\dots j_sj_{s+1}\dots j_r}_{h_1\dots h_sj_{s+1}\dots j_r}=\frac{(n-s)!}{(n-r)!}\delta^{j_1\dots j_s}_{h_1\dots h_s}$$
Where $n$ is the dimension of space.