The curl operator $\vec\nabla\times\mathbb{1}$ can be written as a skew-symmetric 3×3 matrix
$$\mathrm{curl} = \begin{pmatrix}0 & -\partial_z & \partial_y \\ \partial_z & 0 & -\partial_x \\ -\partial_y & \partial_x & 0\end{pmatrix}$$
(in Cartesian coordinates). Since $\mathrm{curl}\,\mathrm{grad}=0$ and $\mathrm{div}\,\mathrm{curl}=0$, the Eigenvector1 to the Eigenvalue 0 is $\vec\nabla$. And by calculating $\det(\mathrm{curl}-\lambda)\stackrel!=0=-\lambda(\lambda^2+\Delta)$ one obtains that the other two Eigenvalues (with opposite signs) both satisfy $\lambda_\pm^2=-\Delta$, i.e. they are the two roots of the negative Laplacian.
Since $\mathrm{curl}^2=\mathrm{grad}\,\mathrm{div}-\Delta$, the Eigenvectors $\vec f_\pm$, satisfying $\mathrm{curl}^2\vec f_\pm = -\Delta\vec f_\pm$, they have a constant divergence. More precisely, since they must be orthogonal to $\vec\nabla$ due to the different Eigenvalues, they can be written as $\vec f_\pm = \vec g_\pm\times\vec\nabla$. Since the $\vec f_\pm$ are also orthogonal to each other, one can obtain
$$\vec f_\pm = \pm\lambda_\pm^{-1}\,\vec f_\mp\times\vec\nabla$$
where the factor $\pm\lambda_\pm^ {-1}$ was chosen for symmetry and consistency. But what is an analytical (non-recursive) expression for them?
You may already have noticed that the $\lambda_\pm$ are Dirac operators, so I wouldn't be too surprised if an answer included spinors and $\partial\!\!/$, even though that is in 4D… In fact, it is most likely that this is the case, since apart from $\vec\nabla$, the vector of $\gamma$ (or Pauli) matrices would not set any preferred direction.
1 In the sense that for any scalar function $s(\vec x)$, $\vec\nabla s$ is an Eigenvector of $\mathrm{curl}$ with Eigenvalue 0, i.e. $\forall s:\mathrm{curl}\vec\nabla s = 0\cdot\vec\nabla s$. It's just like stating $[x_i,\partial_j] = \delta_{ij}$ while this only makes sense when acting on something.
Best Answer
Be $f(x,y,z)=\mathrm e^{\alpha x+\beta y+\gamma z}$ and $v$ an eigenvector of $\begin{pmatrix} 0 & -\gamma & \beta\\ \gamma & 0 & -\alpha\\ -\beta & \alpha & 0\end{pmatrix}$. Then $f(x,y,z)v$ is an eigenvector of $\rm curl$.
Proof: $\partial_x f(x,y,z) = \alpha f(x,y,z)$, $\partial y f(x,y,z) = \beta f(x,y,z)$, $\partial z f(x,y,z) = \gamma f(x,y,z)$. Therefore the "curl matrix" acts on $f(x,y,z)v$ as if all partial derivatives were replaced by the corresponding factor, which gives the matrix above. Thus, we arrive at an ordinary eigenvector equation, and $v$ (and thus each multiple of $v$) is an eigenvector by assumption.
I'm not sure if those are all eigenvectors, though.