Let me give more details. Let $M^{ij} := B^i \wedge B^j$. We can think of $M$ as a matrix with 4-forms as entries. The space of even forms is a commutative algebra, so we can work with $M$ as if it were a real or complex matrix, say. In particular, we can take its trace (which will be a 4-form): $T = \delta_{ij} M^{ij}$, where as in the question the Einstein summation convention is in force. We can then decompose $M$ into a traceless part we shall call $M_0$ and a part containing the trace:$$M^{ij} = M_0^{ij} + \frac{1}{N} T \delta^{ij},$$where I assume that $M$ is an $N\times N$ matrix. If you take the trace of this equation, you find that $M_0$ is indeed traceless. Its explicit form is given by solving that equation for $M_0$, but we do not need it.
Now let $\Phi_{ij}$ be a symmetric traceless matrix. This means that $\delta^{ij} T_{ij} = 0$. Contracting with $M$ we find$$\Phi_{ij} M^{ij} = \Phi_{ij} M_0^{ij}.$$In other words, $\Phi$ never sees the trace of $M$ and hence if you have a lagrange multiplier term in an action functional of the form$$\int \Phi_{ij} M^{ij}$$this is really equal to$$\int \Phi_{ij} M_0^{ij}$$and hence the resulting Euler-Lagrange equation is $M_0^{ij} = 0$.
If $B^i$ are 2-forms, then $B^i \wedge B^j$ is symmetric, not skewsymmetric. Since $\Phi_{ij}$ is traceless, only the traceless part of $B^i \wedge B^j$ that couples to $\Phi$. So I see nothing wrong with the equation you find in the papers.
The reason you take $\Phi$ to be traceless is that the trace is already contained in the second term in the action.