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.