Solving $\nabla \times \mathbf{b} = \mathbf{b} \times \mathbf{a}$
Suppose we are given a fixed vector field $\mathbf{a}$. I am interested in the problem of determining a vector field $\mathbf{b}$ such that $$\nabla \times \mathbf{b} = \mathbf{b} \times \mathbf{a}.$$ This has another interpretation. Suppose $\alpha$ and $\beta$ are the 1-forms dual to $\mathbf{a}$ and $\mathbf{b}$. The above equation can be written as $$d \beta = \beta \wedge \alpha,$$ and so we can interpret this problem as finding, for a fixed 1-form $\alpha$, a foliation $\mathcal{F} = \text{ker}\ \beta$ such that $\alpha$ determines the Godbillon-Vey class of $\mathcal{F}$.
It seems unlikely to me that a solution always exists, but I have been unable to prove anything beyond the obvious fact that we must have $\mathbf{b} \cdot \nabla \times \mathbf{b} = 0$, and that $\mathbf{b} \cdot \nabla \times \mathbf{a} = \mathbf{a} \cdot \nabla \times \mathbf{b}$ (take the divergence), which implies (Mark's comment) that $\mathbf{b}$ is orthogonal to $\nabla \times \mathbf{a}$.
Are there any known results about such equations, or techniques one could use to construct a solution other than crunching through the PDEs for each component?
Here is a small set of solutions. Let $p$ and $q$ be any numbers and $c$ any constant vector. Then $$ b(x) = |x|^p x\times (x\times c), \qquad\qquad a(x) = (p+3)\frac{x}{|x|^2}+qb(x) $$ is a solution. Here $x=(x_1,x_2,x_3)$ and $|x|$ is Euclidean norm.
Edit: your comment about $a$ rules out gradients, but I'll mention anyway more solutions when $a = \nabla g$. Using $\nabla\times(h\nabla f) = \nabla h\times\nabla f$ you get $$ \nabla\times(e^{-g}\nabla f) = e^{-g}\nabla f \times \nabla g. $$