I hope this is an appropriate question for this forum. If not, I apologize. Before stating my question (which may be found at the end of this post), I will attempt to provide sufficient context. 

I am studying eigenfunctions of the Laplacian (both Dirichlet and Neumann cases) on a specified domain $D\subseteq \mathbb{R}^2$ which has piece-wise smooth boundary. Next, consider the billiard map 
$$\phi : M \to M,$$
where $M$ is the subset of the unit cotangent bundle $S^\ast(\partial D)$ containing only the inward pointing directions. Let $p$ be a given function that is invariant with respect to the billiard map. i.e, $$p:M \to \mathbb{C}$$ is such that
$$
p\circ \phi = p \quad \text{on } M.
$$
By quantizing $p$, we obtain a pseudo-differential operator $P_h = p(x, hD)$. I have shown that $P_h$ shares a complete (in $L^2(D)$) collection of eigenfuntions with the Laplacian $\Delta$. In particular, $P_h$ commutes with the Laplacian $\Delta$.

> Is there an intuitive reason why the operator I obtained commutes with the Laplacian? Can anyone provide information regarding the significance (or applications) of such an operator?