Don't use indices. Use a sophisticated and efficient notation for tensor manipulations: through clifford algebra.
Let $\gamma^0, \gamma^1, \gamma^2, \gamma^3$ be basis vectors. Under the clifford product operation, they obey the following:
$$\gamma^\mu \gamma^\nu = \begin{cases}\eta^{\mu \nu} & \mu = \nu \\ -\gamma^\mu \gamma^\nu & \mu \neq \nu\end{cases}$$
The product operation is also associative.
Write $E$ and $B$ in terms of $\gamma^1, \gamma^2, \gamma^3$ as you would write any vectors: $E = E_x \gamma^1 + \ldots$ for instance. Then $F$ can be written
$$F = \gamma^0 E - \gamma^1 \gamma^2 \gamma^3 B$$
Define $\epsilon = -\gamma^0 \gamma^1 \gamma^2 \gamma^3$. Then the tensor equation $\epsilon^{\alpha \beta \gamma \delta} F_{\alpha \beta} F_{\gamma \delta}$ can be written
$$\epsilon^{\alpha \beta \gamma \delta} F_{\alpha \beta} F_{\gamma \delta} = \langle \epsilon FF\rangle_0$$
The brackets denote that we're looking for the scalar part of this expression. This is equivalent to $\epsilon \langle FF \rangle_4 = \epsilon (F \wedge F)$ in the common shorthand.
Now we can use the associativity of the product to our advantage. For brevity, let $i = \gamma^1 \gamma^2 \gamma^3$, such that $\gamma^0 i = \epsilon$, and see that
$$\langle \epsilon FF \rangle_0 = \langle \epsilon (\gamma^0 E - iB)(\gamma^0 E - iB) \rangle_0$$
Multiply out the terms:
$$\langle \epsilon FF \rangle_0 = \langle \epsilon (\gamma^0 E \gamma^0 E -\gamma^0 E iB - iB \gamma^0 E + iBiB )\rangle_0$$
Use the properties of the product to deduce that $E\gamma^0 = -\gamma^0 E$, $\gamma^0 B = -B \gamma^0$, $Ei = iE$ and $iB = Bi$, and $ii = 1$. The result reduces to
$$\langle \epsilon FF \rangle_0 = \langle \epsilon (-EE - 2 \epsilon EB + BB) \rangle_0$$
Finally, deduce that $\epsilon \epsilon = -1$, $EE = -|E|^2$ and $BB = -|B|^2$ to get
$$\langle \epsilon FF \rangle_0 = \langle \epsilon (|E|^2 - |B|^2) + 2 E \cdot B \rangle_0$$
The only term that survives the angled brackets is the $2 E \cdot B$ term. Of course, staring us right in the face is the other well-known invariant, $E^2 - B^2$. That is one advantage of this method: both invariants can be derived from the same approach, and at the same time. Moreover, we never have to break into components to evaluate the expression. We merely use the commutation properties of basis elements in the clifford algebra--an algebra already widely used in quantum mechanics.
For a proof that the only ${\rm SO}(N)$ invariant tensors are products of $\delta_{ab}$'s and Levi-Civita symbols see M. Spivak, A Comprehensive Introduction to Differential Geometry (second edition) Vol. V, pp. 466-481. The number of pages required for the argument shows that it is not trivial.
I've just looked up the third edition of Spivak vol V. What is needed is theorem 35 on page 327. This is the section entitled "A Smattering of Classical Invariant Theory." He writes in terms of scalar invariants, but of course an invariant tensor becomes a scalar invariant when contracted with enough vectors, and any such scalar invariant arises from an invariant tensor.
Best Answer
You can use Young tableaux/diagrams and the permutation group to figure out the symmetries of the general rank-3 tensor. The spaces correspond to the partitions of the rank:
3=3:
One 20 dimensional total symmetric subspace.
3=2+1:
Two 20 dimensional mixed symmetry subspaces.
3=1+1+1:
One 4 dimensional totally antisymmetric subspace:
$$ A_{\alpha\beta\gamma} = \frac 1 6 [T_{\alpha\beta\gamma} + T_{\beta\gamma\alpha} + T_{\gamma\alpha\beta} - T_{\gamma\beta\alpha} - T_{\beta\alpha\gamma} - T_{\alpha\gamma\beta}] $$
That is the only antisymmetric thing you can make according to Schurl-Weyl theory.
To find the dimensions, I used the Hook Length Formula (summing of the boxes $x$ in a diagram $Y(\lambda)$) for the Young diagram corresponding to the integer partition:
$$ {\rm dim}\pi_{\lambda} = \frac {n!}{\prod_{x\in Y}{\rm hook}(x)}$$
If you consider 3 dimensions ($n=3$), you get ${\rm dim} = 1$, that is the standard Levi-Civita symbol $\epsilon_{ijk}$.
If you set $n=4$, the result is ${\rm dim} = 4$.
That means $A_{\alpha\beta\gamma}$ transforms like a 4-vector.
So, the only antisymmetric part of a rank-3 tensor in Minkowski space rotates like a 4-vector, which means it is not invariant and is not a candidate to be Levi-Civita like.
Meanwhile, the dimensions of the 3 other irreducible spaces are all 20--which are certainly not scalars, and thus not candidates to be Levi-Civita like.
Note that if you consider rank-4 tensors, the partitions are as follows:
4=4:
35 dimensional and symmetric.
4=3+1:
Three 45-dimensional mixed symmetry spaces.
4=2+2:
Two 20-dimensional mixed symmetry spaces.
4=2+1+1:
Three 15-dimensional mixed symmetry spaces.
4=1+1+1+1:
One total antisymmetric 1 dimensional space, which is proportional to the Levi-Civita symbol $\epsilon_{\mu\nu\sigma\lambda}$.
In summary, the answer is "No", and the reason why has to do with the representations of the symmetric group on 3-letters. You partition the rank=3, use the Robinson-Schensted correspondence to associate that partition with irreducible representations of the permutation group. (The Young Diagrams make this step a snap). Then, Schur-Weyl duality associates those with irreducible subspaces of and rank-N tensor (signed permutations of indices). Finally, the Hook Length formula tells you the dimensions of those subspaces.
The Levi-Civita symbol needs to be invariant (e.g., dimension 1, like a scalar) and it need to be totally antisymmetric in all indices--and that simply did not exist for rank 3 in 4 dimensions.