This follows from $[X\wedge Y,Z] = X\wedge [Y,Z] \pm [X,Z]\wedge Y $ (see [1])

This does not generalize to $p$-vector fields with $p>2$. One can compute that if $X = X_1 \wedge X_2 \wedge X_3$ then $[X,X] = 0$ always.

This might work. I will think about it.

Thank you for your suggestion. As I answered @WillieWong, requiring $u$ to be a gradient is too restrictive. A generalization of this condition is to require $\delta u = g$, that is, $\partial u / \partial x^i - \partial u/ \partial x_j = f_{ii}$, where $f_{ij} = - f_{ji}$. I could then try to study the operator $P$ augmented with the extra equations. Do you know if there in some sense in which this operator could be considered elliptic?

@WillieWong this is an interesting suggestion, but I think that it restricts too much the classes of solution that I can consider. Locally this is equivalent than considering that the rotational vanishes. The obvious generalization would be to set $δu = g$..

I also though about using the Hodge decomposition and fix the codifferential of $u$ and try to use boundary conditions to determine the harmonic part, but I do not think that this approach is possible if you allow $g$ to depend on $u$

Also, you are right that $f$ is not a top form, but a fiber bundle morphism from $E \to M$ to $\Lambda^k(M) \to M$, so $f(\cdot, u(\cdot))$ is the actual volume form. I only look for locally defined solution, so the manifold is completely irrelevant.

@WillieWong You are right about the EDIT, I will remove it.

What is the deffinition of the Laplace-Beltrami operator? Does your manifold has some Riemannian metric? If so, does it satisfy some compatibility condition with the contact structure (e. g. is M a Sassakian manifold)?

Is there one that has the property that the even numbers have measure 1/2? Or, more generaly, the mesasure of n $\mathbb{N}$ is $1/n$?