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 operatior 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$.
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$, and, expanding both sides, one finds that one must also 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. What is 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$ is a linear total differential equation, 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. Let's 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 not then 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 = \tfrac12\ a'(u)$ and $K = -\tfrac12\ a''(u)$. One also has $L = a(u) -\tfrac12\ a''(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 formula give a solution to the intertwining equation. If $a$ is positive and periodic on the $u$-line with some period, then one gets a solution on the torus.
One can also get solutions on the $2$-sphere: If $a$ is smooth and satisfies $a(u_0) = a(u_1) = 0$ 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, and this gives a family of solutions to the intertwining equation.
To get $\widehat{\Delta^1}$ to be the Bochner Laplacian, of course, you want $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 and be done for the case when you don't make quite as restrictive an assumption about $\widehat{d^\ast}$, but I won't go into that unless someone is interested.