I thought a little bit about your question, which is phrased a little more generally than I like, but I decided to think about it with the restrictions that $\widehat{\Delta^0}$ be a differential operator and $\widehat{d^\ast}$ be $d^\ast$ plus a lower-order (i.e., zeroth order) differential operator. I also decided to think, not of the most general perturbation of $\widehat{\Delta^1}$, but of perturbations of the form $\widehat{\Delta^1} = \nabla^\ast\nabla+L$, where $L$ means scalar multiplication by a function $L$ on $M$.
The first thing to notice is that, since the principal symbol $\sigma^1_\xi$ of $\widehat{\Delta^1}$ is just scalar multiplication by $|\xi|^2$ and the principal symbol of $\widehat{d^\ast}$ is $\alpha\mapsto \xi\cdot\alpha$, it follows from the symbol calculus that the principal symbol of $\widehat{\Delta^0}$ must also be scalar multiplication by $|\xi|^2$, i.e., the differential operator $\widehat{\Delta^0} - \Delta^0$ must be of order at most $1$. However, it must also be self-adjoint, and this implies that it must be of order $0$, i.e., that $\widehat{\Delta^0} = \Delta^0 + H$ (where '$H$' means scalar multiplication by a smooth function $H$ on $M$).
Now, $\widehat{d^\ast}\alpha = d^\ast\alpha + \phi\cdot\alpha$ for some $1$-form $\phi$ on $M$. Under the given assumptions, it is not hard to see that the equation
$$
\widehat{d^\ast}\widehat{\Delta^1} = \widehat{\Delta^0}\widehat{d^\ast}
$$
implies that $\widehat{\Delta^0} = \Delta^0 + H$ for some function $H$ on $M$.
(This is what I just explained above.)
Now, the operator $E = \widehat{d^\ast}\widehat{\Delta^1} - \widehat{\Delta^0}\widehat{d^\ast}$, when expanded out, is of first order (not second order, as you might have expected). Looking at the principal symbol of $E$ and setting this equal to zero gives $4$ equations, and these are equivalent to $\nabla\phi = f\ g$, where, for simplicity, I have set $f = \tfrac12(L-K-H)$. Taking the covariant derivative of both sides of $\nabla\phi = f\ g$ and using the definition of $K$, one gets $df = -K\ \phi$. Substituting this back into the expression for $E$, it reduces to a $0$-th order operator which turns out to be $E(\alpha) = (2f\phi +dK - dL)\cdot\alpha$. Setting this equal to zero gives $d(|\phi|^2 + K - L) = 0$, since $d\bigl(|\phi|^2\bigr) = 2f\phi$ (which is a consequence of $\nabla\phi = f\ g$). Thus, $L = K + |\phi|^2 - c$ for some constant $c$.
Thus, one finds that one must have the identities
$$
\nabla \phi = f\ g\qquad\text{and}\qquad df = -K\ \phi
$$
for some function $f$ on $M$ while
$$
L = K + |\phi|^2 - c\qquad\text{and}\qquad H = |\phi|^2 -2f - c
$$
for some constant $c$.
Obviously, there is always the trivial solution $(\phi,f) = (0,0)$, which gives the well-known intertwining of the Hodge Laplacian on $1$-forms and $0$-forms. More interesting is the case when there are nontrivial solutions $(\phi,f)$ to the above equations. This puts severe restrictions on the metric $g$, but these can be understood.
Let's assume that $M$ is connected and complete. Then the pair of equations $\nabla\phi = f\ g$ and $df = -K\ \phi$ forms a linear total differential system, so if the pair $(\phi,f)$ vanishes anywhere, it vanishes identically. Let's assume that it does not. The equation $\nabla\phi = f\ g$ implies that $d\phi = 0$ and that the vector field $F$ that is $g$-dual to $\ast\phi$ is a Killing field. Thus, $(M,g)$ is, at least locally, a surface of revolution. Any fixed points of $F$ are isolated elliptic points. Write $\phi = du$ for some function $u$ (at the moment locally defined). Note that the critical points of $u$ are nondegenerate and of index 0 or 2. It is then not hard to show that $|\phi|^2 = a(u)$ for some function $a$ and that, in the region where $a>0$, the metric $g$ takes the form
$$
g = \frac{du^2}{a(u)} + a(u)\ d\psi^2
$$
for some (locally defined) function $\psi$. In these coordinates, one has $\phi = du$, $f = \tfrac12a'(u)$, and $K = -\tfrac12a''(u)$. One also has $L = a(u) -\tfrac12a''(u)- c$ and $H = a(u) - a'(u) - c$.
Conversely, if one starts with a function $a$ positive on some domain on the $u$-line, then the above formulae give a solution to the intertwining equation. If $a$ is positive and periodic on the $u$-line, then one obtains a solution on the torus.
One can also obtain solutions on the $2$-sphere: If $a$ is smooth and satisfies $a(u_0) = a(u_1) = 0$ (where $ u_0 < u_1 $) while $a$ is positive between $u_0$ and $u_1$ and satisfies $a'(u_0) = -a'(u_1) = b >0$, then the metric
$$
g = \frac{du^2}{a(u)} + \frac{4a(u)}{b^2}\ d\theta^2
$$
defines a smooth metric on the $2$-sphere with $(u,\theta)$ as 'polar coordinates' (with $\theta$ being $2\pi$-periodic and $u_0\le u\le u_1$, with $u=u_i$ corresponding to the 'poles' of rotation), and this gives a family of solutions to the intertwining equation.
To get $\widehat{\Delta^1}$ to be the Bochner Laplacian, one must have $L = 0$, which is equivalent to $a''(u) = 2\bigl(a(u)-c\bigr)$. This has 'spherical solutions', such as
$$
a(u) = c - e\ \cosh\bigl(\sqrt{2}\ u\bigr),
$$
where the constants $c$ and $e$ satisfy $c > e > 0$. Thus, there is a nontrivial $2$-parameter family of metrics on the $2$-sphere such that the Bochner Laplacian has the desired intertwining property.
A similar calculation can be done for the more general case when $\widehat{d^\ast}$ is allowed to be somewhat more general (but still underdetermined elliptic), but I won't go into that unless someone is interested. I'll just say that, aside from the metrics above, the only metrics that admit a nontrivial solution when $\widehat{d^\ast}$ is an underdetermined elliptic first order differential operator are the metrics of the form
$$
g = \bigl( a\ (x^2{+}y^2) + 2b\ x + c\bigr)\ \bigl(dx^2 + dy^2\bigr)
$$
where $a$, $b$, and $c$ are constants such that the open set $M$ in the $xy$-plane defined by $a\ (x^2{+}y^2) + 2b\ x + c > 0$ is nonempty. Thus, the metrics that allow this are very restricted.
I'm writing about the scalar-valued equation, so the solution is $u \colon \Omega \to \mathbb{R}$ and $\Omega \subset \mathbb{R}^d$, $d \geq 1$. Maybe this is of some use.
Basic lecture notes on $p$-Laplace equation, mathematical aspects: http://www.math.ntnu.no/~lqvist/p-laplace.pdf . See also the bibliography.
The $p$-Laplace equation is a prototype of nonlinear (or quasilinear) elliptic PDE and has many properties that resemble those of the 2-Laplace equation.
Variational justification: The solutions of 2-Laplace equation minimize the energy
$$\int_\Omega |\nabla u|^2 \text{d} x$$
in the space $H^1 (\Omega) = W^{1,2} (\Omega)$ with fixed Dirichlet boundary conditions.
Solutions of the $p$-Laplace equation minimize the energy
$$\int_\Omega |\nabla u|^p \text{d} x$$
in the space $W^{1,p} (\Omega)$ with fixed Dirichlet boundary conditions.
One possible physical interpretation is conductivity of electricity. In your situation there should also be some power-law behaviour.
Recall the Ohm's law, which states that current flux $j$ is proportional to differences in electric potential $\nabla u$ (I assume constant conductivity); $$-j = \nabla u.$$
By Kirchhoff's law you have $\nabla \cdot j = 0$ when there are no sources or sinks of electricity.
Combine these and you have the Laplace equation
$$-\Delta u = 0.$$
The Ohm's law is only an approximation; in reality, you can have complicated non-linear relations there. One possible relation is of power-law type, where
$$-j = |\nabla u|^{p-2}\nabla u,$$
which leads to the $p$-Laplace equation. This power law relation has been observed in some materials near the temperatures where the material becomes superconductive; there $p$ is a function of temperature.
On history: I have a faint memory of someone saying that the origin of $p$-Laplace equation is in (non-linear) fluid dynamics. I have not checked this out. I guess Ladyzhenskaja would be a likely author. Perhaps investigate there?
Best Answer
Disclaimer
The below argument applies directly to the previous version of the question. For the current version the basic argument still hold, using that $$ \| \nabla u\|_2 = \langle u, -\Delta u\rangle $$ and so behaves nicely for eigenfunctions of the Laplacian.
To make the argument also work for $C^\infty_c$ functions simply apply a cut-off near the boundary. This is a small perturbation and so you just need to add some small $\epsilon$ terms to the whole thing.
The added requirement that $\int u = \int v$ is of no consequence since the constructed example below has no zero-frequency component, and that the insertion of $\nabla$ in the new version of the question kills this component anyway.
End disclaimer
With the given information it seems unlikely.
Let $b,g,h$ be three eigenfunctions of $-\Delta$ with eigenvalues $\lambda_b \gg \lambda_g \gg \lambda_h$, normalized so $\|b\|_2 = \|g\|_2 = \|h\|_2 = 1$. Let $\tau_0 = 0$, $\tau_1 = 1/\lambda_g$, and $\tau_2 = 1 / \lambda_h$.
Let $u = a_g g$ and $v = b + a_h h$ where $a_g$ and $a_h$ are scalars. We have that $$ \begin{align} f(\tau_0) &= a_g^2 - (1 + a_h^2) \\ f(\tau_1) &= \frac{a_g^2}{4} - \frac{1}{(1 + \frac{\lambda_b}{\lambda_g})^3} - \frac{a_h^2}{(1 + \frac{\lambda_h}{\lambda_g})^3} > \frac{a_g^2}{4} - \frac{\lambda_g^3}{\lambda_b^3} - a_h^2 \\ f(\tau_2) &= - \frac{a_h^2}{8} - \frac{1}{(1 + \frac{\lambda_b}{\lambda_h})^3} + \frac{a_g^2}{(1 + \frac{\lambda_g}{\lambda_h})^2} < - \frac{a_h^2}{8} + \frac{a_g^2 \lambda_h^2}{\lambda_g^2} \end{align}$$ It is clear that for very large $t$, $f(t) > 0$ (since $t^{-1}$ decays slower than $t^{-3/2}$). So to produce an almost counterexample (only thing missing is that $u$ and $v$ are not in $C^\infty_0(\Omega)$), we just need to arrange
$$\begin{align} a_g^2 &< a_h^2 + 1 \\ a_g^2 & \geq 4 a_h^2 + \frac{4 \lambda_g^3}{\lambda_b^3} \\ a_g^2 & \leq \frac{\lambda_g^2}{8 \lambda_h^2} a_h^2 \end{align} $$
This can be satisfied by, for example, $a_h^2 = 1/6$ and $a_g^2 = 1$ provided the ratio of eigenvalues $24\lambda_g^3 < \lambda_b^3$ and $48 \lambda_h^2 < \lambda_g^2$.