sign error proving product rule for the Laplacian on a product of Riemannian manifolds Given two Riemannian manifolds $M$ and $N$, of dimension $m$ and $n$ respectively, the product manifold $M\times N$ has a Riemannian structure, and there is a Laplacian operator $\Delta_{M\times N}$ on $\Lambda^k(M\times N)\cong \bigoplus_j \Lambda^jM\otimes \Lambda^{k-j}N.$ We expect the Laplacian to obey a product rule identity $\Delta_{M\times N}\omega\wedge\eta=(\Delta_M\omega)\wedge\eta + \omega\wedge(\Delta_N\eta),$ where $\omega\in\Lambda^jM$ and $\eta\in\Lambda^{k-j}N.$
We've been trying to prove this identity using the identities for the exterior derivative and Hodge star operator on products: $d(\omega\wedge\eta)=d\omega\wedge\eta + (-1)^j\omega\wedge d\eta$ and $*(\omega\wedge\eta)=(-1)^{(m-j)(k-j)}*\omega\wedge*\eta$.
So far we've been unable to get the signs to come out right, terms that should cancel aren't, and terms that should combine don't. We're hoping another set of eyeballs will find the sign error. Or perhaps the formula itself is wrong, or one of our assumptions is wrong?
Instead of posting the general proof, which is very lengthy and nearly illegible, I thought I would post instead a computation of a special case in $\mathbb{R}^3$. Let $f(x)g(y,z)\,dy\in\Lambda^1\mathbb{R}^3\cong \Lambda^0\mathbb{R}\otimes\Lambda^1\mathbb{R}^2.$ 


*

*$\Delta(fg\,dy)=(*d*d)(fg\,dy) + (d*d*)(fg\,dy)=(*d*)(g\partial_xf\,dx\wedge dy +f\partial_zg\,dz\wedge dy) + (d*d)(fg\,dz\wedge dx)
=(*d)(g\partial_xf\,dz-f\partial_zg\,dx) + (d*)(f\partial_yg\,dy\wedge dz\wedge dx)=*(g\partial^2_xf\,dx\wedge dz + \partial_xf \partial_yg\,dy\wedge dz-f\partial_y\partial_zg\,dy\wedge dx-f\partial^2_zg\,dz\wedge dx) + d(f\partial_yg)=-g\partial^2_xf\,dy + \partial_xf \partial_yg\,dx+f\partial_y\partial_zg\,dz-f\partial^2_zg\,dy+\partial_xf\partial_yg\,dx+f\partial^2_yg\,dy+f\partial_z\partial_yg\,dz$


Now the Laplacian of a 1-form on $\mathbb{R}^2$ is given by $\Delta(g\,dy)=(\partial^2_y+\partial^2_z)g\,dy$ and the Laplacian on 0-forms on $\mathbb{R}$ is just $\Delta(f)=\partial^2_xf$. So the result I'm expecting is $\Delta(fg\,dy)=(\Delta f)g\,dy+f(\Delta g\,dy)=g\partial^2_xf\,dy+f(\partial^2_y+\partial^2_z)g\,dy.$ If I could flip the sign of the first and fourth terms in 1, and also flip either the second or fifth so they cancel, and also flip either the third or seventh so they cancel, then the result would be proved. 
However we can find no sign error which would justify doing so. What we have instead is $\Delta(fg\,dy)=-(\partial^2_xf)g+f(\partial^2_y-\partial^2_z)g\,dy+(\partial_xf\partial_yg\,dx+\partial_xf\partial_yg\,dx)+(f\partial_y\partial_zg\,dz+f\partial_z\partial_yg\,dz)$ How can I make this calculation work out? 
Maybe some of my Hodge stars are incorrect? I'm putting $*(dx)=dy\wedge dz,*(dy)=dz\wedge dx,*(dz)=dx\wedge dy,*(dx\wedge dy)=dz,*(dy\wedge dz)=dx,*(dz\wedge dx)=dy,*(dx\wedge dy\wedge dz)=1.$ Those formulas are standard on $\mathbb{R}^3$. But could they somehow be different because I'm decomposing $\mathbb{R}^3$ as a product $\mathbb{R}\times \mathbb{R}^2$, and each product factor has its own orientation?
And though I haven't posted it here due to length, we've also worked out the derivation for arbitrary an $k$-form $\omega\wedge\eta$ on $M\times N$, and I have similar problems there; terms that should combine don't, terms that should cancel don't. Can someone help me find some more minus signs in this calculation?
 A: Edit: As mentioned in the comments, the Laplacian on Riemannian manifolds is defined not by $\Delta=*d*d+d*d*$, but $\Delta=\delta d+d\delta$, where the coderivative $\delta$ is defined by $\delta=(-1)^{nk+n+1}*d*$ for a $k$-form in $n$ dimensions. In 3D, this is the same as $\delta=(-1)^k*d*$.
I want to point out an even simpler case which perhaps highlights where you are going wrong.
The laplacian evaluated on $f(x)$ (Note that in the next few calculations the $(-1)^k$ is symbolic and is acting as an operator, i.e. noncommutative, because I am waiting to see what rank the form is before I evaluate it):
$$\Delta f=(-1)^k*d*df + d(-1)^k*d*f=(-1)^k*d*\partial_xf\,dx + d(-1)^k*df\,dx\wedge dy\wedge dz$$
$$=(-1)^k*d\partial_xf\,dy\wedge dz + d(-1)^k*\partial_xf\,dx\wedge dx\wedge dy\wedge dz$$
$$=(-1)^k*\partial^2_xf\,dx\wedge dy\wedge dz=(-1)^0\partial^2_xf=\partial^2_xf$$
Nothing fancy here, just what you'd expect. Now let's try $f(x)\,dy$:
$$\Delta f\,dy=(-1)^k*d*df\,dy + d(-1)^k*d*f\,dy$$
$$=(-1)^k*d*\partial_xf\,dx\wedge dy - d(-1)^k*df\,dx\wedge dz$$
$$=(-1)^k*d\partial_xf\,dz + d(-1)^k*\partial_xf\,dx\wedge dx\wedge dz$$
$$=(-1)^k*\partial^2_xf\,dx\wedge dz=-(-1)^1\partial^2_xf\,dy=\partial^2_xf\,dy$$
So the missing minus sign is coming from this mysterious "codifferential" object. Let's try the full calculation now:
$$\Delta(fg\,dy)=((-1)^k*d*d)(fg\,dy) + (d(-1)^k*d*)(fg\,dy)$$
$$=((-1)^k*d*)(g\partial_xf\,dx\wedge dy +f\partial_zg\,dz\wedge dy) + (d(-1)^k*d)(fg\,dz\wedge dx)$$
$$=((-1)^k*d)(g\partial_xf\,dz-f\partial_zg\,dx) + (d(-1)^k*)(f\partial_yg\,dy\wedge dz\wedge dx)$$
$$=(-1)^1*(g\partial^2_xf\,dx\wedge dz + \partial_xf \partial_yg\,dy\wedge dz-f\partial_y\partial_zg\,dy\wedge dx-f\partial^2_zg\,dz\wedge dx) + d(-1)^0(f\partial_yg)$$
$$=g\partial^2_xf\,dy - \partial_xf \partial_yg\,dx-f\partial_y\partial_zg\,dz+f\partial^2_zg\,dy+\partial_xf\partial_yg\,dx+f\partial^2_yg\,dy+f\partial_z\partial_yg\,dz$$
$$=(\partial^2_xf)g\,dy +f(\partial^2_y+\partial^2_z)g\,dy$$
and all is well.
A: Isn't the Laplacian $\Delta = d\delta + \delta d$?  Wikipedia lists $\delta = (-1)^k *^{-1} d *$.  I think $k$ is meant to be the grade of the form being acted on, but I'm not sure.
Edit: the clumsiness of this minus sign on $\delta$ is, to me, one strike against differential forms in general.  Geometric calculus does a much better job of treating this differential operator on the same footing as the exterior derivative.
In geometric calculus, we just have $\nabla$.  When acting on a multivector field $A$, the exterior derivative is $\nabla \wedge A$, and the codifferential ("interior derivative") is denoted $\nabla \cdot A$.  You can add these together and treat them like separate components of a larger multivector.  We say $\nabla \wedge A + \nabla \cdot A \equiv \nabla A$, and that's that.
Here's an example in 3d:  let's consider the 3d field that you wanted to look at.  This would be $A = f(x) g(y,z) e^y$.  The vector derivative is then
$$\nabla A = g \frac{\partial f}{\partial x} e^x \wedge e^y + f \frac{\partial g}{\partial z} e^z \wedge e^y + f \frac{\partial g}{\partial y}$$
This combines $d$ and $\delta$ into one operator.  A big point of convenience is that the Laplacian is then just $\Delta = \nabla^2$, in all coordinate systems, vector Laplacian or scalar.  And we can compute $\nabla^2 A$ with no issue.
$$\begin{align*}\nabla^2 A &= \nabla (\nabla A) \\&= g \frac{\partial^2 f}{\partial x^2} e^y - \frac{\partial g}{\partial y} \frac{\partial f}{\partial x} e^x + f \frac{\partial^2 g}{\partial z^2} e^y - f \frac{\partial^2 g}{\partial z \partial y} e^z \\& \quad + \frac{\partial f}{\partial x} \frac{\partial g}{\partial y} e^x + f \frac{\partial^2 g}{\partial y^2} e^y + f \frac{\partial^2 g}{\partial y \partial z} e^y \\ &= g \frac{\partial^2 f}{\partial x^2} e^y + f \left(\frac{\partial^2 g}{\partial y^2} + \frac{\partial^2 g}{\partial z^2} \right) e^y \\ &= g\nabla^2 f + f \nabla^2 g\end{align*}$$
Exactly as you expected.  All the minus signs are automatically taken care of by $e^x \wedge e^y = -e^y \wedge e^x$ and contractions.
As with differential forms, the equality of mixed partials gives us identities about closed forms.  $\nabla \wedge \nabla \wedge F = 0$ for any $F$, and $\nabla \cdot (\nabla \cdot F) = 0$ for the same reason (these are equivalent to $d^2 = \delta^2 = 0$).  This makes the calculation of the above much simpler.
It is my hope that seeing what geometric calculus does well (treating the coderivative) will give some insight into how painful it is to treat in differential forms.  Still, differential forms is the standard by which any other formalism that tries to capture multivariate calculus is judged by.  I hope this merely convinces you that the coderivative is, well, difficult and clumsy.
Finally, I tried my hand at proving the overall goal through GC.  I imagine it's as simple as saying that $\nabla_{M+N} = \nabla_M + \nabla_N$ and then multiplying.  Doing that gives
$$\nabla^2_{M+N} = \nabla_M^2 + \nabla_N^2 + \nabla_M \nabla_N + \nabla_N \nabla_M$$
The symmetric product will naturally kill any bivector terms.  What concerns me is that there will be mixed scalar cross terms if the spaces $M$ and $N$ are not mutually orthogonal.  Perhaps this is implicit or obvious.  Still, the question of non-orthogonal coordinate systems is what concerns me here.
