Pfaffian of a skew symmetric $n\times n$ matrix $A$ with $n=2m$ even is defined as
$$Pf(A)=\frac{1}{2^m m!} \sum_{i_1,...,i_n} \epsilon^{i_1...i_n} a_{i_1i_2}\cdot \cdot \cdot a_{i_{n-1}i_n}, $$
where $\epsilon^{i_1...i_n}$ are Levi Civita symbols:
$$\epsilon^{i_1...i_n}:= \begin{cases}
+1&\text{if } (i_1,...,i_n) \text{ is an even permutation of } (1,...,n) \\
-1&\text{if } (i_1,...,i_n) \text{ is an odd permutation of } (1,...,n)\\
0& \text{else }
\end{cases}$$
As direct consequence Leibniz formula for determinant becomes
$$ det(B)=\sum_{i_1,...,i_n} \epsilon^{i_1...i_n}b_{1i_1}\cdot \cdot \cdot b_{ni_n},$$
multiplied both sides by $\epsilon^{j_1...j_n}$ we have
$$ \epsilon^{j_1...j_n}det(B)=\sum_{i_1,...,i_n} \epsilon^{i_1...i_n} b_{j_1i_1}b_{j_2i_2} \cdot \cdot \cdot b_{j_{n-1}i_{n-1}} b_{j_ni_n}. $$
Now recall that $(B^tAB)_{\textbf{ij}}=\sum_{k=1}^n \sum_{l=1}^n b_{k\textbf{i}} a_{kl} b_{l\textbf{j}}$ and notice that we are free to choose the name for indices $k,l$ in the sum. They will become $j_{s-1},j_{s}$ in equation below.
\begin{align}
2^m m!Pf(B^tAB)&=\sum_{i_1,...,i_n} \epsilon^{i_1...i_n} (B^tAB)_{i_1i_2}\cdot \cdot \cdot(B^tAB)_{i_{n-1}i_n}\\
&=\sum_{i_1,...,i_n} \epsilon^{i_1...i_n} \sum_{j_1,...,j_n}\big(b_{j_1i_1}a_{j_1j_2}b_{j_2i_2}\big) \cdot \cdot \cdot \big(b_{j_{n-1}i_{n-1}}a_{j_{n-1}j_n}b_{j_ni_n}\big)\\
&=\sum_{j_1,...,j_n} \sum_{i_1,...,i_n} \big[ \epsilon^{i_1...i_n} b_{j_1i_1}b_{j_2i_2} \cdot \cdot \cdot b_{j_{n-1}i_{n-1}} b_{j_ni_n}\big] (a_{j_1j_2} \cdot \cdot \cdot a_{j_{n-1}j_n})\\
&=\sum_{j_1,...,j_n} \epsilon^{j_1...j_n}det(B) (a_{j_1j_2} \cdot \cdot \cdot a_{j_{n-1}j_n})=det(B) 2^m m! Pf(A)
\end{align}
This proof is taken from Michael Spivaks Comprehensive Introduction to Differential Geometry. vol 5-Publish or Perish (1999).