I am trying to understand parametrices for the solution operator $G_t = \sin(t\sqrt{\Delta})/\sqrt{\Delta}$ to the wave equation $$(\partial_{tt} + \Delta)u=0, ~~~~~~~ u_0 =0, ~~~~~~\partial_tu_0 = f$$ on a manifold with smooth boundary, i.e. if $W(t, x, y)$ is the Schwartz kernel of the operator $W_t$, there should be distributions $R_j(t, x, y)$ supported in the "light cone" $\{t^2\geq d(x, y)^2\}$ such that the difference $$W(t, x, y) - \sum_{j=0}^N R_j(t, x, y)$$ is of regularity $C^k$.
Much work has been done by Taylor and Melrose who constructed the distributions $R_j$ as Fourier-Airy integrals. However, the assumption is always that the boundary have either only positive curvature or only negative curvature (glancing versus gliding rays). Also the analysis is not as implicit as in the boundary-less case, e.g. as in Bär-Gauduchon-Pfäffle, who express the $R_j$ very explicit and geometrically in terms of Riesz-Distributions.
On the other hand, the papers by Taylor and Melrose I found are more than 20 years old, and I wonder if there has been any progress in this direction.
Have people constructed a parametrix making no restriction on the boundary? Is there a more geometric approach?
