Tensorial Proof that $\det(F^{\mu\nu})=(\vec{E}\cdot\vec{B})^2$? I am trying to understand why 
$$\det(F^{\mu\nu})=(\vec{E}\cdot\vec{B})^2\tag{1}.$$ 
Of course one can just calculate the determinant of $F^{\mu\nu}$ expressed as a matrix with components given in terms of $E_x, E_y, ...$, etc., but I am looking for something a bit more insightful. In particular, I understand that $\vec{E}\cdot\vec{B}$ can be written as 
$$\vec{E}\cdot\vec{B}=-\frac{1}{8}\epsilon_{\mu\nu\rho\sigma}F^{\mu\nu}F^{\rho\sigma}.\tag{2}$$
This seems promising because then
$$(\vec{E}\cdot\vec{B})^2=\frac{1}{64}\epsilon_{\mu\nu\rho\sigma}\epsilon_{\alpha\beta\gamma\delta}F^{\mu\nu}F^{\rho\sigma}F^{\alpha\beta}F^{\gamma\delta},\tag{3}$$
which looks awfully like the expression for the determinant of $F$ in terms of Levi-Civita tensors,
$$\det(F^{\mu\nu})=\frac{1}{4!}\epsilon_{\mu\nu\rho\sigma}\epsilon_{\alpha\beta\gamma\delta}F^{\mu\alpha}F^{\nu\beta}F^{\rho\gamma}F^{\sigma\delta}.\tag{4}$$
But I can't figure out how to connect these two expressions. Since $F$ is antisymmetric, one can flip the two indices of any single $F$ at the cost of a minus sign, but I need a way to permute indices between different $F$'s. I'm also not sure how a factor of 3 could possibly enter in to change the 1/64 into a 1/24. 
 A: The two expressions are only superficially similar; they have to be by dimensional analysis. The index structure is completely different, and I don't think you can convert one to the other without essentially undoing everything you did to get to the first result, then inserting the standard proof of the second result. Such a proof would be unenlightening since it'd just be unnecessarily complicated.
So let me just give a simple derivation of the result you want, directly. The determinant is
$$\det F^{\mu\nu} = \epsilon_{\mu\nu\rho\sigma} F^{\mu 0} F^{\nu 1} F^{\rho 2} F^{\sigma 3}$$
There are naively $4!$ terms, but the only nonzero terms are those where $\mu \neq 0$, $\nu \neq 1$, $\rho \neq 2$, $\sigma \neq 3$. That is, we need to count the number of derangements of a 4-element set. It is easy to show by casework there are $9$.
Each of these $9$ terms is quadratic in $\mathbf{E}$, since there are two indices equal to zero, and hence quadratic in $\mathbf{B}$. Moreover, the sum of all $9$ terms is a tensorial invariant. There are only two independent invariants: $E^2 - B^2$ and $\mathbf{E} \cdot \mathbf{B}$. We can't use the first one, because otherwise our terms wouldn't all have the same degree in $E$ and $B$. Then the answer must be proportional to $(\mathbf{E} \cdot \mathbf{B})^2$, but since $(\mathbf{E} \cdot \mathbf{B})^2$ has $9$ terms, they must simply be equal.
You might complain I used components, but I had to because your expression is not tensorial. The fields $\mathbf{E}$ and $\mathbf{B}$ are not Lorentz tensors, but rather a way of writing components of $F_{\mu\nu}$. You can't expect to prove a statement about components without expanding in components.
A: OP's eq. (2) is (minus) the Pfaffian ${\rm Pf}(F)$, and OP's eq. (4) is the determinant ${\rm Det}(F)$. The sought-for relation (1) follows because the square of the Pfaffian is the determinant: 
$$ {\rm Pf}(F)^2~=~{\rm Det}(F). \tag{A}$$
A proof of eq. (A) is given in my Math.SE answer here.
A: The determinant on the left is invariant under (proper) Lorentz transformations, so you can easily calculate it in a frame of reference where $E^2=E^3=0$, it equals $(E^1 H^1)^2=(\vec{E}\vec{H})^2 $. As $\vec{E}\vec{H} $ is an invariant under proper Lorentz transformations, this proves your formula in a general case.
A: Quoting myself from an old PhysicsForums post

If I am not mistaken, (someone check my math)
  you can calculate the principal invariants of a tensor:
  $$F^a{}_a,\quad F^a{}_{[a} F^b{}_{b]},\quad F^a{}_{[a} F^b{}_{b}F^c{}_{c]},\quad F^a{}_{[a} F^b{}_{b} F^c{}_{c}F^d{}_{d]}$$
  in 4-dimensions. 
  You'd get the trace (sum of the eigenvalues), sum of products-of-pairs of eigenvalues, sum of products-of-triples, and finally, the product-of-the-4-eigenvalues (the determinant).
  
  For a real antisymmetric matrix, the eigenvalues are imaginary.
  So, only those with sums of products of even-numbers-of-eigenvalues will be nonzero.

(Higher combinations like $F^a{}_{[a} F^b{}_{b} F^c{}_{c}F^d{}_{d}F^e{}_{e}F^f{}_{f]}$ vanish since it would have more than 4 factors being antisymmetrized). So, all but two are identically zero.
At some point in your calculation, you would probably have to make use of the epsilon-delta identities [sometimes written as determinants with delta-entries] (like this contracted combination $\epsilon_{abmn}\epsilon^{cdmn}=-4\delta_a{}^{[c}\delta_b{}^{d]}$ (up to sign conventions, based on signature and dimensionality)). Together with the antisymmetry of $F_{ab}=F_{[ab]}$, I think the index-gynmastics will work out.
