Why are we interested in operators that share a basis of eigenfunctions? 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 bounded 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. The billiard map models the behaviour of a free particle in the space $D$. Suppose we are base point $x\in \partial D$ and a unit direction vector $w$ pointing inwards. Then a free particle starting at x and travelling in direction $w$ will eventually hit the boundary $\partial D$. Let $x^\prime$ be the point of incidence and suppose that $w^\prime$ is the new direction of the particle upon hitting the boundary. Then  $\phi(x, w) = (x^\prime, w^\prime)$. The Billiard map can also be described in terms of the billiard flow. That is, suppose that $\varphi_t : S^\ast D \to S^\ast D$ is the billiard flow. The billiard flow $\varphi_t$ solves the equation
$$
\partial_t\varphi_t(x, \omega) = \begin{pmatrix}
0 & 1\\
0 & 0
\end{pmatrix}\varphi_t(x, w).
$$
Then
$$
\phi(x, w) = \varphi_{\tau(x,w)}(x, w)
$$
where
$$
\tau(x, w) = \inf\{t > 0: \varphi_t(x, w) \in M\}.
$$
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.
$$
Supposed that the function $p$ is defined on all of $\mathbb{R}^2\times \mathbb{R}^2$. By quantizing $p$, we obtain a pseudo-differential operator $P_h = p(x, hD)$. That is,
$$
P_hu(x) = \frac{1}{\left(2\pi h\right)^2}\int_{\mathbb{R}^2}\int_{\mathbb{R}^2} e^{\frac{i}{h}\langle x-y, \xi\rangle} p(x, \xi) u(y)\,\mathrm{d}{y}\mathrm{d}{\xi}
$$
for every sufficiently nice function $u$.
I have shown that $P_h$ shares a complete (in $L^2(D)$) collection of eigenfuntions with the Laplacian $\Delta$. We also note that, $P_h$ commutes with the Laplacian $\Delta$.

Is there an intuitive reason why the operator I obtained shares an orthonormal basis of (Dirichlet/Neumann) eigenfunctions with the Laplacian? Can anyone provide information regarding the significance (or applications) of such an operator? 

 A: I'm not sure what precisely you are looking for but let me try the guess an answer on the soft end (basic intuition stuff) of what you are asking. I'm not qualified for the hard end of this question, for this you want an expert on partial differential operators and microlocal analysis. I'm in the applied/physics corner.
To me the relationship of billiard constructions to PDEs is most evident when one considers a dynamical case, e.g. the d'Alembertian over the domain. This then is the physical case of the dynamics of a drum-head (for a compact domain). One way to attach this problem was to hope that separation of spacial (Laplacian) and temporal (second partial derivative with respect to the time signature) and then a further separation in over the domain with respect to some coordinate system yield something we can say a lot about. The general strategy here is "global" (in the sense of the size of the domain) Fourier decomposition of the correct eigenfunctions defined by the solutions of separated ODEs. For a set of simple domains (rectangle, circle, ellipse) this was possible but in general it's not. 
The reason why the d'Alembertian is nice here for intuition is because it more strongly motivates a more local approach. Rather than study the behavior "globally" over the range of the domain, what if we model local behavior and observe "evolution" of the solution. Early attempts at this on hyperbolic PDEs is the theory of characteristics. If you inject a local non-trivial initial condition (say a impulsive distribution), the solution will propagate in a certain way (along characteristics). Think of the initial impulse position to be an infinite set of infinitessimal billiard balls, and think of the d'Alembertian imposing the dynamics all these balls have to follow, and the initial directions of each ball corresponds to the characteristic cone. Then you will get a billiard decomposition of this local disturbance. If you just study the trajectories of these billiard balls you get the billiard maps with their reflection laws. If however you also study either local value of displacement, you end up studying what is known the fundamental solution of the PDE (the impulse response of that PDE). The billiard maps in this setting captures the directions of the characteristics for one initial location. The full solution would be recovered by convolution with an arbitrary initial function over the domain. Billiards at least intuitively are dynamical objects which is why I like to think of them in the hyperbolic case like this.
Given that we are in a linear operator theory here we have some nice transformation and decomposability properties. So we can reduce to the Laplacian by separating out the time component. In the classical PDE theory this is separating temporal oscillations from oscillatory shapes (eigenvalues of complex exponentials over time and eigenfunctions over the domain shape). So the study of the Laplacian here is merely the reduction to only study the eigen"shapes" of say a drum-head. So when we use billiard-type decompositions we are in a sense only looking for the stationary solutions and not the temporal dynamics. But we can still do this kind of decomposition regardless. Though we consider billiard maps to be a flow the purpose is to represent the Laplacian alone.
The goal of course is for your the full collection of functions over billiard maps to contain all the information that is the solution of the laplacian over the domain. So the eigenfunctions (in the re-composition) have to be identical!
But the well-known cases show that one has to be careful here. The typical eigenfunctions of a single billiard map is the complex exponential. Yet the eigenfunction of the Laplacian over a circular domain is Bessel. So your precise treatment of eigenfunctions in the re-composition of all billiard maps with all the local behavior has to be equal to the Bessel function to solve the problem. But this is an easy case because we have the Bessel ODE and lots (though not all) useful properties so we can compare. In many non-trivial cases all we can do is move forward and make sure all our steps are correct. (So when you say that you have the eigenfunction of the Laplacian over your domain, I'm a bit confused how you got it, but there is lots of detail missing). It turns out that local singularities of the billiard map bundle relate to "phase" contributions. We count the number of these with Keller-Maslov (or more generalized) indices, and the nature of the phase change is subject to local integration around the singularity.
Basic topic areas and important contributors (apologies for any omissions!) for related literature:
Physics: Geometric optics, wave optics, semi-classical physics.


*

*Keller's work on geometric optics and correction to quantization conditions. Maslov's work. Ballian Bloch. Gutzwiller. See also Berry on the relation of geometric and wave optics and Nye's beautiful book on the topic.


Mathematics:


*

*Microlocal analysis. Index theories. Singularity theory, particularly in the context of Lagrangian and Legendrian manifolds. I also recommend some Arnold for a highly geometric view of what happens at singularities (his work on wavefronts and caustics). Integration theory of singularities in particular Picard-Lefschetz.


Video:
Zeldich recently gave a recently talk on his result that mild elliptic perturbations of the circle are spectrally determined (not isospectral) contains a lot of related background:
https://www.youtube.com/watch?v=sL73RYikETw
