Determinant of an $8 \times 8$ antisymmetric matrix

I want to show that the square root of the determinant of the following antisymmetric matrix

is given by

I saw this equation in the paper Supplemental Material of Four-Dimensional Quantum Hall Effect with Ultracold Atoms, but it didn't provide a specific deduction.

I have tried Mathematica and Matlab to calculate the determinant directly. But as for Mathematica, the processing time is too long to wait. And though Matlab could give the result very fast but it couldn't show them all because the result is too long. Thus, I guess the best way to solve this problem might be mathematical deduction rather than computer. After all the matrix is antisymmetric, and it's quite special in many aspects.

Second, let us denote by $M$ the matrix in question and write $M=J+R$, where $R$ is $M$ with all $\pm 1$'s replaced by zeros. By linearity of the Pfaffian, we can write $$\mathrm{Pf}\,M=\sum_S \epsilon(S,S')\mathrm{Pf}\,J_{S'}\mathrm{Pf}\,R_S,$$ where the sum is over the set of all even-sized subsets of $I:=\{1,\dots ,8\}$, $S'=I\setminus S$, and $\epsilon(S,S')$ is the signature of the permutation in which elements of $S'$ follow the elements of $S$ in the increasing order. In fact, $\mathrm{Pf}\,J_{S}=\pm 1$ if each set $\{1;5\}$,$\{2;6\}$,$\{3;7\}$,$\{4;8\}$ is contained in either $S$ or $I\setminus S$, and zero otherwise.
Now, $S=\emptyset$ contributes $1$, the first term in the answer, and $S=I$ contributes $\mathrm{Pf}\,R=\mathrm{Pf}\,B\cdot \mathrm{Pf}\,\tilde{\Omega}$, which is the last term. If $S$ is a four-element set (e. g., $\{1;5;2;6\}$), then $R_S$ is block-diagonal with two blocks of size $2$, and the Pfaffian is just the corresponding product (in the above case, $-B_{12}\tilde{\Omega}^{56}$). Summing over all these cases gives the second term in your formula. If $S$ is two-element or six-element, then $R_S$ will be block-diagonal with blocks of odd size, and hence the Pfaffian will be zero.
You are looking at the Pfaffian $Pf_A$ of the matrix $$A=\begin{pmatrix} B & -I_4 \\ I_4 & \Omega \end{pmatrix}.$$ It is a polynomial in the entries, of total degree $4$. By Schur complement formula, we have $$\det A=\det B\det(\Omega+B^{-1})=\det(I_4+B\Omega).$$ Let us introduce a parameter $t$ and form $A(t)$ with the matrices $tB$ and $\Omega$ instead. Then $Pf_{A(t)}^2=\det(I_4+tB\Omega)$ and $Pf_{A(t)}$ is a quadratic polynomial in $t$. Denoting $\lambda_j$ the eigenvalues of $B\Omega$, we have $$\det(I_4+tB\Omega)=1+t\sum_j\lambda_j+t^2\sum_{j<k}\lambda_j\lambda_k+O(t^3)$$ from which we infer $$Pf_{A(t)}=1+\frac12t\sum_j\lambda_j+\frac12t^2\sum_{j<k}\lambda_j\lambda_k-\frac18t^2(\sum_j\lambda_j)^2.$$ We thus have $$Pf_A=1+\frac12{\rm Tr}(B\Omega)+\frac12\sum_{j<k}\lambda_j\lambda_k-\frac18(\sum_j\lambda_j)^2,$$ where the last two sums form a polynomial in $B$ and $\Omega$.
I now claim that the multiplicities of the eigenvalues of $B\Omega$ are even. This because the characteristic polynomial satisfies $$\chi(X)=\det(B\Omega-XI_4)=\det B\det(\Omega-XB^{-1})=(Pf_B\cdot Pf_{\Omega-XB^{-1}})^2$$ is the square of a polynomial. Let me call $\lambda,\mu$ the eigenvalues, which are therefore double. We write the last sums $$\frac12\sum_{j<k}\lambda_j\lambda_k-\frac18(\sum_j\lambda_j)^2=\lambda\mu=(\det B\Omega)^{\frac12}=(\det B\det\Omega)^{\frac12}=\pm Pf_BPf_\Omega.$$
We therefore obtain $$Pf_A=1+\frac12{\rm Tr}(B\Omega)\pm Pf_BPf_\Omega.$$ The sign $\pm$ is determined by taking a specific example. For instance if $\Omega=B$ and the eigenvalues of $B$ are $\pm ia,\,\pm ib$, we find easily $$Pf_A=1-a^2-b^2+a^2b^2=1+\frac12{\rm Tr}(B^2)+(Pf_B)^2.$$ Therefore the sign $\pm$ is $+$.