Can the curl operator be generalized to non-3D?
Solution 1:
div, grad, and curl are secretly just the three exterior derivatives in $\mathbb{R}^3$. Said another way, they are the three nontrivial differentials in the de Rham complex
$$0 \to \Omega^0(\mathbb{R}^3) \xrightarrow{d_0} \Omega^1(\mathbb{R}^3) \xrightarrow{d_1} \Omega^2(\mathbb{R}^3) \xrightarrow{d_2} \Omega^3(\mathbb{R}^3) \to 0$$
of $\mathbb{R}^3$. In particular, grad is secretly $d_0$, curl is secretly $d_1$, and div is secretly $d_2$. I say "secretly" because there is some additional funny business going on involving the Hodge star and the musical isomorphisms.
So in $n$ dimensions there are $n$ generalizations of div, grad, and curl $d_0, ... d_{n-1}$ which satisfy $d_{i+1} \circ d_i = 0$ (the defining equation of a chain complex).
Solution 2:
In geometric algebra and calculus, the cross product and curl are generalized by the wedge product and exterior derivative.
Let's start with the cross product. Let $u = a \times b$. We usually say $u$ is orthogonal to both $a$ and $b$, or normal to the plane spanned by $a$ and $b$. Geometric algebra works with that plane directly, as a "2-vector" or bivector that we call $a \wedge b$. This is something that can be in any number of dimensions, unlike the cross product.
We compute the exterior derivative of a field $A$ as $\nabla \wedge A$ in the same manner. Together, the exterior derivative with the interior derivative $\nabla \cdot A$ constitute the full, geometric derivative $\nabla A$. This operator has all of the information needed to reconstruct the original $A$ through integration.
Thus, if one knows $\nabla A$ and $\nabla \cdot A$, then one can always compute the curl through $\nabla \wedge A = \nabla A - \nabla \cdot A$.
The relationship to the Laplacian: the exterior derivative and interior derivative obey the following relationships in flat space.
$$\nabla \wedge \nabla \wedge A = \nabla \cdot (\nabla \cdot A) = 0$$
This follows as a consequence of the equality of mixed partial derivatives. This means that the Laplacian can always be written as
$$\nabla (\nabla A) = \nabla^2 A = \nabla \cdot (\nabla \wedge A) + \nabla \wedge (\nabla \cdot A)$$
This is true for both the scalar and vector Laplacian, although when $A$ is a scalar field, $\nabla \cdot A = 0$ always. The dot should be interpreted as "contracting" or "grade-lowering", and clearly you can't reduce a scalar to something lower-dimensional.
Edit: to build off of what Qiaochu Yuan is saying, let's look at various fields in 3d to understand how this theory captures the usual vector calculus derivatives. Let $E$ be a vector field and $\varphi$ a scalar field. We usually have the following derivatives:
Divergence: $\nabla \cdot E$
Gradient: $\nabla \varphi$
Curl: $\nabla \times E$
The relationship between curl and the wedge derivative is through Hodge duality. This is denoted by multiplication with a pseudoscalar, called $i$. $i$ converts vectors to 2-vectors ("pseudovectors") and vice versa. The canonical relationship between the curl and this derivative is
$$\nabla \times E = -i \nabla \wedge E$$
$\nabla \wedge E$ is what differential forms people call the curl.
What's interesting to note is when the derivative acts on a 2-vector instead. Let such an object be $B$:
$$\nabla \cdot B = \text{vector}, \quad \nabla \wedge B = \text{pseudoscalar (0d vector space)}$$
This latter is what differential forms people tend to call the divergence. To me, this is deceptive: $\nabla \cdot E$ is just as equally the divergence.
Edit edit: you mention complex numbers. Let's talk about them and spinors. You can use this definition of the vector derivative $\nabla$ and its associated operations to work on spinors also. In GA, spinors are just linear combinations of scalars and 2-vectors. Let $\psi(z)$ be such a "spinor" field, or a function of a so-called complex variable. If $\psi = u + iv$ is holomorphic, then it obeys
$$\nabla \psi = \nabla u + (\nabla v) i = 0$$
This may seem somewhat backwards from complex analysis, but this is because $\nabla$ can be identified with $\partial/\partial \bar z$. The benefit here is, of course, that $\nabla$ is valid in all dimensions. You shouldn't be afraid of spinors. They're very useful, and we can treat derivatives of them in the same way we treat fields of other kinds. The spinors of 3d are quaternions, and their usefulness for 3d rotations is well-known.
Finally, let's calculate some derivatives of some sample fields. Let $\varphi$ be a scalar field, $E$ be a vector field, $B$ a bivector field, and $\gamma$ a pseudoscalar field. Let $e_x, e_y, e_z$ be the unit cartesian vectors. Let $e_{zx} = e_z \wedge e_z$, and so on. Let $e^x$ be the unit cotangent vector perpendicular to the $yz$ plane. Note that for Cartesian coordinates, this is equal to $e_x$, but vector derivatives always depend on these one-forms, so this is more general. A quick summary:
$$\newcommand{curl}{\operatorname{curl}}\newcommand{grad}{\operatorname{grad}}\renewcommand{div}{\operatorname{div}}\begin{align*} \varphi &= \varphi(x, y, z) \\ E &= E^x(x, y, z) e_x + E^y e_y + E^z e_z \\ B &= B^{xy} e_{xy} + B^{yz} e_{yz} + B^{zx} e_{zx} \\ \gamma &= \gamma^{xyz}(x, y, z) e_{xyz} \end{align*}$$
Interior derivatives: $$\begin{align*} \nabla \cdot \varphi &= 0 \\ \nabla \cdot E &= \partial_x E^x + \partial_y E^y + \partial_z E^z \\ \nabla \cdot B &= (\partial_z B^{zx} - \partial_y B^{xy}) e_x + (\partial_x B^{xy} - \partial_z B^{yz}) e_y + (\partial_y B^{yz} - \partial_x B^{zx}) e_z \\ \nabla \cdot \gamma &= \partial_x \gamma^{xyz} e_{yz} + \partial_y \gamma^{xyz} e_{zx} + \partial_z \gamma^{xyz} e_{xy}\end{align*}$$
Exterior derivatives:
$$\begin{align*} \nabla \wedge \gamma &= 0 \\ \nabla \wedge B &=( \partial_x B^{yz} + \partial_y B^{zx}+ \partial_z B^{xy}) e_{xyz} \\ \nabla \wedge E &= (\partial_z E^y - \partial_y E^z) e_{yz} + (\partial_x E^z - \partial_z E^x) e_{zx} + (\partial_y E^x - \partial_x E^y) e_{xy} \\ \nabla \wedge \varphi &= \partial_x \varphi e^x + \partial_y \varphi e^y + \partial_z \varphi e^z\end{align*}$$
Edit$^3$: let's talk about Hodge duality. The unit pseudoscalar $i$ turns dots to wedge and wedges to dots. This is true outside of the calculus, too, in a generalized sense: let $b$ be a vector. Then $b \cdot b = b \wedge (bi) i^{-1}$. This is how dot products are often calculated in differential forms (I personally find it unduly circuitous). Applying this to the calculus yields various relations. Let $\varphi, E, B, \gamma$ be defined as before.
$\nabla \cdot \gamma = [\nabla \wedge (\gamma i)] i^{-1}$: this converts the interior derivative of a pseudoscalar field to a gradient of a scalar field, which is then dualized back to a pseudoscalar. This is one way traditional vector calculus converts back and forth between scalars and pseudoscalars, avoiding the difference between the two.
$\nabla \wedge B = [\nabla \cdot (Bi)] i^{-1}$: this is traditionally how differential forms treads the divergence, by considering the exterior derivative of a 2-form.
Edit$^4$: it may be I breezed right past what is of interest here. Let me translate a few things to make apparent a broader point.
$$\begin{align*} \text{Vector calculus:} &\curl \grad \psi = 0 \\ \text{Geometric calculus:} & \nabla \wedge \nabla \wedge \psi = 0 \end{align*}$$
Here, $\nabla \wedge \psi$ is the gradient, and the second $\nabla \wedge$ on the left plays the role of curl.
$$\begin{align*} \text{Vector calculus:} &\div \curl E = 0 \\ \text{Geometric calculus:} & \nabla \wedge \nabla \wedge E = 0 \end{align*}$$
On the other hand, here the divergence and curl both take on the meaning of $\nabla \wedge$. This is why differential forms likes to call them all aspects of $\nabla \wedge$. But both of these identities for divergence, curl, and gradient have exactly the same expression in both differential forms and geometric calculus.
Again, remembering that $\nabla \wedge \nabla \wedge X = 0$ for any kind of $X$ is what allows us to write
$$\begin{align*} \text{Vector calculus:} &\curl \curl E = \grad \div E - \Delta E \\ \text{Geometric calculus:} & -\nabla \cdot (\nabla \wedge E) = \nabla \wedge (\nabla \cdot E) - \nabla^2 E\end{align*}$$
This is an important point: the second applied curl ends up taking the role of $\nabla \cdot$ on the 2-form $\nabla \wedge E$. This is why, ultimately, I think your search for a generalized curl specifically may be somewhat misguided, for curl, divergence, and even gradient take on wildly different meanings depending on what they're acting on, while $\nabla \wedge$ and $\nabla \cdot$ have fixed meanings in geometric calculus as raising (for $\nabla \wedge$) or lowering ($\nabla \cdot$) the dimensionality of the objects they differentiate. It's for this reason that, in geometric calculus, $\nabla \cdot$ is sometimes referred to as the divergence and $\nabla \wedge$ is referred to as the curl, but this lazy (yet useful) terminology doesn't always coincide with vector calculus. Again, the reason that the meanings of divergence, curl, and gradient fluctuate in vector calculus is because of the use of Hodge duality to simplify all objects to scalars and vectors.
Solution 3:
Based on Qiaochu Yuan's answer and Branimir Ćaćić's comments to it, I propose the following generalizations for $n\ge2$ dimensions:
$$\newcommand{curl}{\operatorname{curl}}\newcommand{grad}{\operatorname{grad}}\renewcommand{div}{\operatorname{div}}\begin{array}{rl} \grad &:= \sharp\circ d_0,\\ \div &:= \ast\circ d_{n-1}\circ\ast\circ\flat,\\ \curl^2 &:= \sharp\circ\ast\circ d_{n-2}\circ\ast\circ d_1\circ\flat \end{array}$$
using the path $$\Gamma(\mathbb R^n)\xrightarrow{\flat}\Omega^1(\mathbb R^n)\xrightarrow{d_1}\Omega^2(\mathbb R^n)\xrightarrow{\ast}\Omega^{n-2}(\mathbb R^n)\xrightarrow{d_{n-2}}\Omega^{n-1}(\mathbb R^n)\xrightarrow{\ast}\Omega^1(\mathbb R^n)\xrightarrow{\sharp}\Gamma(\mathbb R^n)$$
I am not yet sure how to express $\Delta$ and therefore can't check whether $\curl^2=\grad\div-\Delta$. Since $d_{i+1}\circ d_i = 0$, one obtains $\curl^2\grad=0$ and $\div\curl^2=0$, which is necessary (but not sufficient) for the criteria $\curl\grad=0$ and $\div\curl=0$ from the question. These require
$$\curl = \sharp\circ\ast\circ d_{n-2}\circ A\circ d_1\circ\flat$$ for some $A$ which, by comparing $\curl^2$, must obey $$A\circ d_1\circ\ast\circ d_{n-2}\circ A = \ast$$ or $$d_1\circ\ast\circ d_{n-2} = A^{-1}\circ\ast\circ A^{-1}.$$ So, can $A$ be $d_1^{-1}$ and $d_{n-2}^{-1}$ at the same time for $n\neq3$?