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.
Best Answer
Q: Does anybody know how to numerically overcome this pseudospectral effect?
The key idea is "normal ordering". Rewrite the problem in such a way that annihilation operators $a$ appear to the right of creation operators $a^\ast$. In this particular case, first notice that $H$ has chiral symmetry, if $\lambda$ is an eigenvalue then also $-\lambda$ is an eigenvalue. We can thus reconstruct the full spectrum of $H$ without sign ambiguities from the spectrum of $H^2$, $$H^2=\begin{pmatrix} 1+a^\ast a&0\\ 0&a^\ast a \end{pmatrix}.$$ We can now safely truncate each block to $N$ states. The spectrum of $H_N^2$ contains 0 as eigenvalue with multiplicity 1.