Wedge product and its relation to determinant. 
Let $\omega_1=\sum_{i=1}^2\sum_{j=1}^2a_{ij}dx_idx_j$,
$\omega_2=\sum_{i=1}^4\sum_{j=1}^4a_{ij}dx_idx_j$,
$\omega_3=\sum_{i=1}^N\sum_{j=1}^Na_{ij}dx_idx_j$.
In each case, I'd like to find $\omega\wedge\omega$ and find the relationship between $\omega\wedge\omega$ and the matrix with entries $a_{ij}$.

I think $\omega_1\wedge\omega_1=0$ because each term will will look like $dx_1dx_2dx_1dx_2$, and each of these is zero because $dx\wedge dx=0$.
For $\omega_2\wedge\omega_2$, the only nonzero terms are permutations of $dx_1dx_2dx_3dx_4$. There are twelve of these, and after rearranging, you get something like $fdx_1dx_2dx_3dx_4$, but what is an easy way to see what $f$ is and its relation to the matrix with entries $a_{ij}$?
and $\omega_n\wedge\omega_n$?
 A: Take
$$
\omega=\sum_{j,k=1}^Na_{jk}{\rm d}x_j\wedge{\rm d}x_k
$$
for instance. Note that $j$ and $k$ are merely indices for summation. Thus by exchanging all $j$'s and $k$'s, we have
$$
\omega=\sum_{j,k=1}^Na_{jk}{\rm d}x_j\wedge{\rm d}x_k=\sum_{j,k=1}^Na_{kj}{\rm d}x_k\wedge{\rm d}x_j.
$$
Yet due to the asymmetry of the wedge product,
$$
{\rm d}x_k\wedge{\rm d}x_j=-{\rm d}x_j\wedge{\rm d}x_k
$$
holds unconditionally. Therefore, the expression goes that
$$
\omega=\sum_{j,k=1}^Na_{jk}{\rm d}x_j\wedge{\rm d}x_k=\sum_{j,k=1}^Na_{kj}{\rm d}x_k\wedge{\rm d}x_j=\sum_{j,k=1}^N\left(-a_{kj}\right){\rm d}x_j\wedge{\rm d}x_k.
$$
Due to the equality of the second and the last term, we have
\begin{align}
\omega&=\sum_{j,k=1}^Na_{jk}{\rm d}x_j\wedge{\rm d}x_k=\sum_{j,k=1}^N\left(-a_{kj}\right){\rm d}x_j\wedge{\rm d}x_k\\
&=\frac{1}{2}\left(\sum_{j,k=1}^Na_{jk}{\rm d}x_j\wedge{\rm d}x_k+\sum_{j,k=1}^N\left(-a_{kj}\right){\rm d}x_j\wedge{\rm d}x_k\right)\\
&=\sum_{j,k=1}^N\frac{a_{jk}-a_{kj}}{2}{\rm d}x_j\wedge{\rm d}x_k.
\end{align}
Therefore, given a $2$-form $\omega$, we are always able to play the above trick, making the coefficient in front of ${\rm d}x_j\wedge{\rm d}x_k$ asymmetric. In other words, given
$$
\omega=\sum_{j,k=1}^Na_{jk}{\rm d}x_j\wedge{\rm d}x_k,
$$
we can always write
$$
\omega=\sum_{j,k=1}^Nb_{jk}{\rm d}x_j\wedge{\rm d}x_k,
$$
where
$$
b_{jk}=\frac{a_{jk}-a_{kj}}{2}=-b_{kj}.
$$
Now, apply the above result to our target:
\begin{align}
\omega\wedge\omega&=\left(\sum_{i,j=1}^Na_{ij}{\rm d}x_i\wedge{\rm d}x_j\right)\wedge\left(\sum_{k,l=1}^Na_{kl}{\rm d}x_k\wedge{\rm d}x_l\right)\\
&=\left(\sum_{i,j=1}^Nb_{ij}{\rm d}x_i\wedge{\rm d}x_j\right)\wedge\left(\sum_{k,l=1}^Nb_{kl}{\rm d}x_k\wedge{\rm d}x_l\right)\\
&=\sum_{i,j,k,l=1}^Nb_{ij}b_{kl}{\rm d}x_i\wedge{\rm d}x_j\wedge{\rm d}x_k\wedge{\rm d}x_l.
\end{align}
Play the index-exchanging trick again:
\begin{align}
\omega\wedge\omega&=\sum_{i,j,k,l=1}^Nb_{ij}b_{kl}{\rm d}x_i\wedge{\rm d}x_j\wedge{\rm d}x_k\wedge{\rm d}x_l\\
&=\sum_{i,j,k,l=1}^Nb_{li}b_{jk}{\rm d}x_l\wedge{\rm d}x_i\wedge{\rm d}x_j\wedge{\rm d}x_k\\
&=\sum_{i,j,k,l=1}^Nb_{kl}b_{ij}{\rm d}x_k\wedge{\rm d}x_l\wedge{\rm d}x_i\wedge{\rm d}x_j\\
&=\sum_{i,j,k,l=1}^Nb_{jk}b_{li}{\rm d}x_j\wedge{\rm d}x_k\wedge{\rm d}x_l\wedge{\rm d}x_i.
\end{align}
Now, reorder the indices as $\left(i,j,k,l\right)$ by using the asymmetry property:
\begin{align}
\omega\wedge\omega&=\sum_{i,j,k,l=1}^Nb_{ij}b_{kl}{\rm d}x_i\wedge{\rm d}x_j\wedge{\rm d}x_k\wedge{\rm d}x_l\\
&=\sum_{i,j,k,l=1}^N\left(-b_{li}b_{jk}\right){\rm d}x_i\wedge{\rm d}x_j\wedge{\rm d}x_k\wedge{\rm d}x_l\\
&=\sum_{i,j,k,l=1}^Nb_{kl}b_{ij}{\rm d}x_i\wedge{\rm d}x_j\wedge{\rm d}x_k\wedge{\rm d}x_l\\
&=\sum_{i,j,k,l=1}^N\left(-b_{jk}b_{li}\right){\rm d}x_i\wedge{\rm d}x_j\wedge{\rm d}x_k\wedge{\rm d}x_l.
\end{align}
This result eventually yields
\begin{align}
\omega\wedge\omega&=\sum_{i,j,k,l=1}^N\frac{b_{ij}b_{kl}-b_{li}b_{jk}+b_{kl}b_{ij}-b_{jk}b_{li}}{4}{\rm d}x_i\wedge{\rm d}x_j\wedge{\rm d}x_k\wedge{\rm d}x_l\\
&=\sum_{i,j,k,l=1}^N\frac{b_{ij}b_{kl}-b_{li}b_{jk}}{2}{\rm d}x_i\wedge{\rm d}x_j\wedge{\rm d}x_k\wedge{\rm d}x_l.
\end{align}
Since $\omega\wedge\omega$ is a $4$-form, in two- and three-dimensional cases, the result must be $0$. From four-dimensions onwards, the result is not necessarily zero, e.g.,
$$
\omega=\alpha{\rm d}x\wedge{\rm d}y+\beta{\rm d}z\wedge{\rm d}w\Longrightarrow\omega\wedge\omega=2\alpha\beta{\rm d}x\wedge{\rm d}y\wedge{\rm d}z\wedge{\rm d}w.
$$
A: In general, given a skew-symmetric matrix $A = (A_{ij})_{i,j=1}^{2k}$, i.e. $A_{ij} = -A_{ji}$, one can build the two-form $\omega_A = \frac{1}{2}\sum_{i,j=1}^{2k} A_{ij}\mathrm dx^i\wedge\mathrm dx^j = \sum_{1\le i < j\le 2k} A_{ij}\mathrm dx^i\mathrm dx^j$. Then $\omega_A^k = \operatorname{Pf}(A)\mathrm dx^1\wedge\cdots\wedge\mathrm dx^{2k}$ is a multiple of the volume form on $\mathbb R^{2k}$, and this multiple is called the Pfaffian $\operatorname{Pf}(A)$. Explicitly, it is calculated by summing over all ways of pairing up the elements of $\{1,2,\cdots,2k\}$, with each pairing $(i_1,j_1),\dots,(1_k,j_k)$ contributing $\pm\prod_{n=1}^k A_{i_nj_n}$, where the sign is that of the permutation sending $\{1,\cdots,2k\}$ to $\{i_1,j_1,\cdots,i_k,j_k\}$. It turns out that $\operatorname{Pf}(A)^2 = \operatorname{det}(A)$; you can check this by expanding both sides, but the easier proof is to note that both sides are independent of the basis used to express $A$, and there is always a basis such that $A$ is block skewsymmetric, i.e. of the form 
$$
\begin{pmatrix}
0 & -\lambda_1 & & & \\
\lambda_1 & 0 & & & \\
& & \ddots & & \\
& & & 0 & -\lambda_k \\
& & & \lambda_k & 0\end{pmatrix}
$$
which has Pfaffian $\lambda_1\cdots\lambda_k$ and determinant $(\lambda_1\cdots\lambda_k)^2$.
Now your question is about $\omega_A\wedge\omega_A$, which will in general be some $4$-form. If you want to find the coefficient of $\mathrm dx^a\wedge\mathrm dx^b\wedge\mathrm dx^c\wedge\mathrm dx^d$, you might as well restrict to the subspace spanned by these four basis vectors. The previous discussion shows that the restriction of $\omega_A\wedge\omega_A$ to this subspace is the Pfaffian of the $4\times 4$-minor of $A$ whose rows and columns are indexed by $\{a,b,c,d\}$. This Pfaffian is a square root of the determinant of this minor. Thus in general the coefficients of $\omega_A\wedge\omega_A$ will be square roots of the determinants of these minors.
If $k=1$, there is no such minor, so $\omega_A\wedge\omega_A = 0$. If $k=2$, there is exactly one such minor, so $\omega_A\wedge\omega_A = \pm\sqrt{\det A}\mathrm dx^1\wedge\cdots\wedge\mathrm dx^4$.
