text:
\section{Introduction}\label{Introduction section}
\noindent The purpose of this paper is to develop a Hadamard-Riesz type parametrix for the wave propagator in bounded planar domains with smooth, strictly convex boundary. We then use this parametrix to produce asymptotic expansions for the distributional wave trace near isolated lengths in the length spectrum. Let $\Omega$ be such a domain and denote by $\Delta$ the Dirichlet Laplacian on $\Omega$. The even wave propagator $E(t)$ is defined to be the solution operator for the wave equation
\begin{align}\label{wave equation}
\begin{cases}
(\d_t^2 - \Delta) E = 0 & (x \in \Omega)\\
E(0) = \text{Id} & (x \in \Omega)\\
\d_t E \big|_{t = 0} = 0 & (x \in \Omega),
\end{cases}
\end{align}
with Dirichlet boundary conditions. In spectral theoretic terms, we can write $E(t) = \cos t \sqrt{-\Delta}$, which is the even part of the half wave propagator $e^{i t \sqrt{-\Delta}}$. In \cite{Ch76}, it is shown that there exist Lagrangian distributions $E_j(t,x,y)$ such that microlocally away from the tangential rays, the Schwartz kernel of $E(t)$ is given by
\begin{align}\label{Chazarain cosine}
E(t,x,y) = \sum_{j = - \infty}^{+\infty} E_j(t,x,y) + C^\infty(\mathbb{R} \times \Omega \times \Omega),
\end{align}
with $E_j$ corresponding to a wave of $j$ reflections at the boudary. The sum in \eqref{Chazarain cosine} is locally finite in time. If we restrict our attention to waves which make $j$ reflections and travel approximately once around the boundary, we have the following explicit parametrix for $E_j(t,x,y)$.
\begin{theo}\label{HRP}
Let $\Omega \subset \mathbb{R}^2$ be a bounded domain with smooth, strictly convex boundary. Then there exists $j_0 = j_0(\Omega)\in \mathbb{N}$ sufficiently large such that the following holds: for all $j \geq j_0$, there exists a tubular neighborhood $U_j$ of the diagonal of the boundary $\Delta \d \Omega \subset \Omega \times \Omega$ such that for all $(x, y) \in U_j$ and $t$ less than but sufficiently close to $|\d \Omega|$,
$$
E_j(t,x,y) = \sum_{\pm} \sum_{k = 1}^8 (-1)^j e^{\pm i \pi / 4} \int_{0}^\infty e^{\pm i\tau(t - \Psi_j^k(x,y))} b_{j, k, \pm} (\tau, x, y) d\tau + C^\infty(\mathbb{R} \times U_j),
$$
microlocally near geodesic loops of rotation number $1/j$. Here, $b_{j, k, \pm} \in S_{\text{cl}}^{1/2}(U_j \times \mathbb{R}^1)$ are classical elliptic symbols of order $1/2$ and the functions $\Psi_j^k(x,y)$ ($1 \leq k \leq 8$) are lengths of the $8$ billiard orbits connecting $x$ to $y$ in $j$ reflections and approximately one rotation (see Theorem \ref{8 orbit lemma}). The principal term in the asymptotic expansion for $b_{j,k,\pm}$ is given in boundary normal coordinates $x = (\mu, \phi)$, $y = (\nu, \theta)$ (see Section \ref{Oscillatory integral representation}) by
\begin{align*}
\frac{\pm |\tau|^{1/2}}{2} f^{1/2}(\mu, \phi) \Bigg|\frac{\d \Psi_j^k}{\d \mu} \frac{\d \Psi_j^k}{\d \theta}\frac{\d^2 \Psi_j^k}{\d \nu \d \phi} &+ \frac{\d \Psi_j^k}{\d \phi}\frac{\d \Psi_j^k}{\d \nu}\frac{\d^2 \Psi_j^k}{\d \theta \d \mu}\\ &-     \frac{\d \Psi_j^k}{\d \mu}\frac{\d \Psi_j^k}{\d \nu}\frac{\d^2 \Psi_j^k}{\d \theta \d \phi}       - \frac{\d \Psi_j^k}{\d \phi} \frac{\d \Psi_j^k}{\d \theta} \frac{\d^2 \Psi_j^k}{\d \nu \d \mu} \Bigg|^{1/2},
\end{align*}
where $f(\mu, \phi)$ is a positive function which is smooth up to the boundary $(\mu = 0)$ and satisfies $f(0, \phi) = 1$.
\end{theo}
\noindent Theorem \ref{HRP} bears a remarkable resemblance to Hadamard's parametrix for the wave propagator on boundaryless manifolds (see \cite{Hadamardparametrix1}, \cite{Hadamardparametrix1}), where the phase functions $\pm \tau(t - \Psi_j^k(x,y))$ are replaced by $\pm \tau(t - r(x,y))$, with $r(x,y)$ the geodesic distance from $x$ to $y$. In that case, for $x,y$ near the diagonal, there is only one geodesic connecting $x$ to $y$. In our setting, Theorem \ref{8 orbit lemma} shows that there are exactly $8$ orbits connecting $x$ to $y$ in $j$ reflections and approximately one rotation.
\begin{rema}
For $j \geq j_0(\Omega)$ and $x, y \in \d \Omega$ sufficiently close to the diagonal, there exist only two orbits connecting $x$ to $y$ in $j$ reflections and approximately one rotation. One is in the clockwise direction and the other is in the counterclockwise direction. In particular, when $x = y \in \d \Omega$, there exists a unique geodesic loop based at $x$ of rotation number $1/j$. Restricting the functions $\Psi_j^k(x,y)$ appearing in Theorem \ref{HRP} to the diagonal of the boundary yield the $j$\textit{-loop function}, which we denote by $\Psi_j(q,q)$.
\end{rema}
\noindent As in the case of boundaryless manifolds, we can use the explicit parametrix in Theorem \ref{HRP} to prove trace formulas. It is known that $E(t)$ has a well defined distributional trace
\begin{align}\label{Wave trace}
\text{Tr} \cos t \sqrt{-\Delta} = \sum_{j = 0}^\infty \cos t \lambda_j,
\end{align}
where $\lambda_j^2$ are the Dirichlet eigenvalues of $-\Delta$. The sum in \eqref{Wave trace} converges in the sense of tempered distributions and has singular support contained in the length spectrum
$$
\text{LSP}(\Omega) = \overline{ \{\text{lengths of periodic billiard trajectories} \} } \cup \{0\},
$$
together with $- \text{LSP}(\Omega)$ (see Section \ref{Computing the wave trace}). Each periodic billiard orbit in $\Omega$ can be classified according to its winding number $m$ and the number of reflections $n$ made at the boundary. Denote the collection of periodic orbits of this type by $\Gamma(m,n)$, normalized so that $m \leq n/2$. $\Gamma(m,n)$ is never empty by a theorem of Birkhoff \cite{Birkhoff}. The length spectrum can be decomposed accordingly as
\begin{align}
\text{LSP}(\Omega) = \bigcup_{m, n \in \mathbb{N}} \text{length}(\Gamma(m,n)) \cup \mathbb{N} |\d \Omega|,
\end{align}
where $|\d \Omega|$ is the length of the boundary. Using the parametrix in Theorem \ref{HRP}, we can prove the following theorem.
\begin{theo}\label{Main theorem}
Assume the conditions from Theorem \ref{HRP} hold and $\Omega$ satisfies the noncoincidence condition:
\begin{align}\label{NCC}
\begin{split}
\text{there exists}\,\,\epsilon_0 > 0\,\, \text{such that} \,\,\bigcup_{\substack{m \geq 2\\ n \geq 1}} \text{length}(\Gamma(m,n)) \cap (|\d \Omega| - \epsilon_0, |\d \Omega|) = \emptyset.
\end{split}
\end{align}
For $j \geq j_0$, define $t_j = \inf_{q \in \d \Omega} \Psi_j(q,q)$ and $T_j = \sup_{q \in \d \Omega} \Psi_j(q,q)$. Then on any sufficiently small neighborhood of $[t_j,T_j]$, $\text{Tr} E(t)$ has the asymptotic expansion
\begin{align*}
(-1)^j \Re \left\{e^{i \pi / 4} \int_{\d \Omega} \int_{0}^\infty e^{i \xi (t- \Psi_j(q,q) )} a^j(q,\xi)  d\xi dq\right\} + C^\infty(\mathbb{R}),
\end{align*}
where $a^{j}(q,\xi)$ is a classical elliptic symbol of order $1/2$ with principal part given by
\begin{align*}
a_0^j(q,\xi) =   {4}\xi^{1/2} \sin \omega_{j, 1}(q,q) \sin^{1/2} \omega_{j,2}(q,q) \left|\frac{\d \omega_{j,1}}{\d q'} (q,q) \right|^{1/2} X(q) \cdot N(q).
\end{align*}
Here, $X(q)$ is the position vector to a boundary point $q$ with respect to a fixed origin and $N(q)$ is the outward unit normal at $q$. The angles $\omega_{j,1}, \omega_{j,2}$ are the initial and final angles respectively of the unique billiard orbit $\gamma_j(q,q')$ which connects nearby boundary points $q$ and $q'$ in $j$ reflections and approximately one counterclockwise rotation. The function $\Psi_j(q,q')$ is the length of $\gamma_j(q,q')$ and its restriction to the diagonal is the $j$-loop function.
\end{theo}
\begin{rema}
As the position vector $X$ is chosen with respect to an arbitrary interior point $p$, we can integrate out this symmetry over any measurable subset of the interior in the variable $p$ to obtain a more invariant formula. In particular, we can integrate over open sets, curves and by a limiting argument, the boundary itself in order to obtain a smooth density on $\d \Omega$.  The noncoincidence condition \eqref{NCC} can be weakened and is in particular satisfied for ellipses (see \cite{GuMe79a}) and nearly circular domains (see \cite{HeZe19}).
\end{rema}
\noindent Theorem \ref{Main theorem} provides an explicit formula for the principal term in the parametrix developed in \cite{MaMe82}. In contrast to the methods employed in \cite{MaMe82}, the proof developed in the remainder of this paper uses Hadamard's variational formula for the wave trace and the explicit parametrix for the wave propagator appearing in Theorem \ref{HRP}. These techniques are of independent interest. Theorem \ref{Main theorem} also provides clarity on a discrepancy in the literature regarding the order of the wave trace (cf. \cite{Popov1994}, \cite{MaMe82}). Note that in Theorem \ref{Main theorem}, no assumptions are made on the nondegeneracy of orbits. If the length spectrum has high multiplicity, the wave trace is in general quite complicated. However, when periodic orbits come in a one parameter family corresponding to a caustic, we have the following:
\begin{theo}\label{Rational caustic theorem}
Suppose $\Omega$, $j \geq j_0$ are as in Theorem \ref{Main theorem} and $\mathcal{C}$ is a caustic for $\Omega$ of rotation number $1/j$. If periodic orbits tangent to $\mathcal{C}$ have length $L_j$, then near $t = L_j$, $\text{Tr} E(t)$ has the leading asymptotic
\begin{align*}
c_j \Re \left\{e^{i\pi/4} (t- L_j - i 0)^{-3/2}\right\},
\end{align*}
where $c_j$ is a wave invariant given by the formula
\begin{align*}
c_j = (-1)^{j+1} {4} \int_{\d \Omega} \sin^{3/2} \omega_{j,1}(q,q) \left|\frac{\d \omega_{j,1}}{\d q'}(q,q)\right|^{1/2} X(q) \cdot N(q) dq.
\end{align*}
In this case, $\omega_{j,1}(q,q) = \omega_{j,2}(q,q)$ is the measure of the angle of incidence for the unique periodic orbit of rotation number $1/j$ based at $q \in \d \Omega$.
\end{theo}
\noindent While KAM theory provides the existence of irrational caustics, it is shown in \cite{KaloshinZhangRationalCaustics} that for a fixed $j$, the set of all smooth convex domains possessing a rational caustic of rotation number $1/j$ is polynomially dense in the variable $1/j$ in the collection of all smooth domains, equipped with the $C^\infty$ topology. Exponential density is also proven in the analytic category. As ellipses satisfy the noncoincidence condition \eqref{NCC} (see \cite{GuMe79a}) and are known to be completely integrable with confocal conic sections as caustics, we have the following corollary to Theorem \ref{Rational caustic theorem}.
\begin{coro}\label{Ellipse corollary}
For an ellipse $\Omega$ given by
$$
\Omega = \left\{(x,y): \frac{x^2}{a^2} + \frac{y^2}{b^2} \leq 1\right\},
$$
and $L_j \in \text{LSP}(\Omega)$ sufficiently close to $|\d \Omega|$ corresponding to the length of billiard orbits of rotation number $1/j$, the wave invariants in Theorem \ref{Rational caustic theorem} are given by
\begin{align*}
c_j = \int_0^{2\pi} \frac{(-1)^{j+1} 2 ab \sin \omega_j \sqrt{a^2 \cos^2 \phi + b^2 \sin^2 \phi} d\phi }{\sqrt{\cos \omega_j (a^2 \sin^2 \phi + b^2 \cos^2\phi) (b^2 + (a^2 - b^2) \sin^2 \phi) G(\zeta_j) \sqrt{1 - k_{\zeta_j}^2 \sin \phi}}}.
\end{align*}
Here, $\zeta_j \in [0,b)$ is the parameter of the confocal ellipse
$$
\mathcal{C}_{\zeta_j} = \left\{(x,y): \frac{x^2}{a^2 - \zeta_j^2} + \frac{y^2}{b^2 - \zeta_j^2} = 1\right\},
$$
to which the orbits of length $L_j$ are tangent and $k_{\zeta_j}$ is given by
$$
k_{\zeta_j}^2 = \frac{a^2 - b^2}{a^2 - {\zeta_j}^2}.
$$
The parameter $\zeta_j$ depends on the rotation number ${1/j}$ and is defined implicitly by the equation
\begin{align*}
\frac{1}{j} = \frac{F(\arcsin \zeta_j/b; k_{\zeta_j})}{2 K(\zeta_j)},
\end{align*}
where $F(s; k)$ is the elliptic integral
\begin{align*}
\int_0^{s} \frac{d \tau}{\sqrt{1 - k^2 \sin^2 \tau }}
\end{align*}
and $K(\zeta_j) = F(\pi/2; k_{\zeta_j})$. The function $G(\zeta_j)$ is defined by
\begin{align*}
\frac{-k_{\zeta_j}^2}{(a^2 - \zeta_j^2)} \int_0^{2\pi} \frac{\sin^2\tau d\tau}{(1 - k_{\zeta_j}^2 \sin^2 \tau)^{3/2}} + (2j+2) \frac{d}{d \zeta^2} F(\arcsin \zeta / b; k_{\zeta})\Big|_{\zeta = \zeta_j},
\end{align*}
and $\omega_j = \omega_j(\phi)$ are the angles of reflection for orbits tangent to $\mathcal{C}_{\zeta_j}$, given implicitly by the equation
\begin{align*}
\zeta_j^2 = \sin^2 \omega_j (b^2 + (a^2 - b^2)\sin^2 \phi).
\end{align*}
\end{coro}

\begin{rema}
Analagous formulas to those appearing in Theorem \ref{Main theorem}, Theorem \ref{Rational caustic theorem} and Corollary \ref{Ellipse corollary} can also be proved for the Neumann and Robin wave traces. The formulas are less succinct but can be easily reproduced by following the proof of Theorem \ref{Main theorem} and combining it with the Hadamard type variational formulas for Robin boundary conditions in the author's previous work \cite{Vig18}.
\end{rema}

\subsection{Schematic Outline}
The proofs of Theorem \ref{Main theorem}, Theorem \ref{Rational caustic theorem} and Corollary \ref{Ellipse corollary} use techniques from \cite{Vig18}, which are reviewed throughout the paper. In Section \ref{Background}, we review relevant background on the inverse spectral problem, i.e. determining geometric information from the Laplace spectrum. The first step in the proof of Theorem \ref{Main theorem} is to generalize the Hadamard-Riesz type parametrix for the wave propagator constructed in \cite{Vig18} for ellipses to arbitrary bounded domains with strictly convex, smooth boundary as in Theorem \ref{HRP}. This requires a dynamical classification (Theorem \ref{8 orbit lemma}) of the cardinality and structure of all billiard orbits connecting interior points with a fixed number of reflections, analagous to Lemma 5.2 in \cite{Vig18}. Section \ref{Billiards} introduces language from dynamical systems necessary to describe the billiard (or broken bicharacteristic) flow, which later appears in the canonical relations of the wave propagator. The proof of Theorem \ref{8 orbit lemma} is relegated to Section \ref{Proof of 8 orbit lemma}, where it is broken up into several intermediate lemmas. This material is largely independent from the rest of the paper and is of separate interest from the perspective of dynamical billiards, irrespective of applications to spectral theory. In Section \ref{A parametrix for the wave propagator}, the length functionals corresponding to these orbits allow us to cook up explicit phase functions which parametrize the canonical relations for the wave propagator $e^{i t \sqrt{-\Delta}}$, which is a Fourier integral operator microlocally away from the tangential rays. This parametrix is microlocalized near approximate geodesic loops, which are particularly useful when integrating the wave kernel over the diagonal in the trace formulas which follow. Construction of an explicit parametrix for the wave propagator in the interior appears to be new in the literature, with the exception of the author's prior work \cite{Vig18} in the special case of an ellipse. In Section \ref{Computing the wave trace}, we then use this parametrix to obtain a singularity expansion for a variation of the wave trace corresponding to radial dilation about an interior point, following the Hadamard type variational formulas in \cite{HeZe12} and \cite{Vig18}. The variational formulas are of interest because they allow us to restrict our study of the wave kernels to their behavior at the boundary. This technique for studying the wave trace at a fixed domain by variational methods also appears to be new in the literature. In Section \ref{Principle symbol computation}, the honest principal symbol of the wave trace is computed in a specific coordinate system by integrating its variation in the radial variable. This is essentially the fundamental theorem of calculus combined with some rescaling symmetries of the wave equation. In Section \ref{Integrating the principle symbol}, it is shown that radially integrating subprincipal terms coming from the variational wave trace doesn't contribute to the principal symbol. As the order of the principal symbol computed in Section \ref{Principle symbol computation} appears to contradict other works in the literature, Section \ref{Check on the order} provides an auxilliary confirmation via stationary phase that the order acheived in this paper is indeed correct. To our knowledge, this paper contains the first explicit formulas for the principal symbol of the wave trace associated to a convex billiard table.

\section{Background}\label{Background}
\noindent The inverse spectral problem has a long history, dating back to Kac in 1966 \cite{Kac66}, who asked the famous question ``can one hear the shape of a drum?'' Mathematically, this corresponds to uniquely determining a domain $\Omega$ from the spectrum of its Dirichlet, Neumann or Robin Laplacian. For bounded domains, the spectrum is purely discrete, consisting of eigenvalues $\lambda_j^2$ satisfying
\begin{align}\label{inverse problem}
\begin{cases}
-\Delta u_j = \lambda_j^2 u_j & x \in \Omega\\
Bu_j = 0 & x \in \d \Omega,
\end{cases}
\end{align}
where $u_j$ are smooth eigenfunctions on $\Omega$ and $B$ is either the restriction operator (Dirichlet boundary conditions), normal differentiation (Neumann boundary conditions) or normal differentiation plus a prescribed function on $\d \Omega$ (Robin boundary conditions). A variety of approaches in local and global harmonic analysis have been taken to prove partial results in the direction of \cite{Kac66}. One particularly useful strategy is to use the wave group to deduce spectral information about the underlying geometric space, usually a Riemannian manifold. The motivation behind this approach stems from Duistermaat and H\"ormander's propagation of singularities theorem, which says that singularities of solutions to the wave equation propagate along (possibly broken) bicharacteristics, which are lifts to $T^* (\mathbb{R} \times \Omega)$ of geodesic or billiard orbits. As linear waves can be superimposed, constructive interference is most pronounced along geodesics which are traversed infinitely often, i.e. periodic orbits.  On the trace side, this is reflected in the Poisson relation:
\begin{align}\label{Poisson relation}
\text{SingSupp Tr} e^{i t \sqrt{-\Delta} } \subset \pm \text{LSP} (\Omega) \cup \{0\},
\end{align}
where the lefthand side is the distributional trace of the half wave propagator (see Section \ref{Computing the wave trace}) and the righthand side is the length spectrum of $\Omega$ (the closure of all lengths of periodic geodesic or billiard orbits together with $\{0\}$). This is due to \cite{AndersonMelrose} in the case of bounded planar domains. In particular, the formula \eqref{Poisson relation} generalizes the Poisson summation formula on the torus $\mathbb{R}^n/\mathbb{Z}^n$ from elementary Fourier analysis (see \cite{Ur98}, for example).\\
\\
Asymptotic formulas near the singularities are given by the Selberg trace formula in the case of hyperbolic surfaces \cite{Sel56}, the Duistermaat-Guillemin trace theorem for general smooth manifolds under a dynamical nondegeneracy condition \cite{DuGu75}, and a Poisson summation formula for strictly convex bounded planar domains due to Guillemin and Melrose \cite{GuMe79b}. However, since these trace formulas involve sums over all periodic orbits of a given length, it is theoretically possible that the contributions of distinct orbits having the same length could cancel out and the wave trace is actually smooth near a point in the length spectrum. We say that the length $L \in \mathbb{R}$ of a periodic orbit $\gamma$ is simple if up to time reversal ($t \mapsto -t$), $\gamma$ is the unique periodic orbit of length $L$. Without length spectral simplicity, there is no way to deduce Laplace spectral information from the length spectrum alone. It is shown in \cite{PeSt92} that generically, smooth convex domains have simple length spectrum associated to only nondegenerate periodic orbits. In that case, the following theorem holds:
\begin{theo}[\cite{GuMe79b}, \cite{PeSt17}]
Assume $\gamma$ is a nondegenerate periodic billiard orbit in a bounded, strictly convex domain with smooth boundary and $\gamma$ has length $L$ which is simple. Then near $L$, the even wave trace has an asymptotic expansion
\begin{align}
\text{Tr} \cos t \sqrt{-\Delta} \sim \Re \left\{ a_\gamma (t - L + i0)^{-1} + \sum_{k = 0}^\infty a_{\gamma k} (t - L + i0)^{k} \log(t - L + i0) \right\},
\end{align}
where the coefficients $a_{\gamma k}$ are the wave invariants associated to $\gamma$.
\end{theo}
\noindent Calculations of the wave invariants associated to dynamically convenient orbits have proved extremely useful in the inverse spectral problem associated to \eqref{inverse problem}. For example, the case of rotationally symmetric metrics on $S^2$ is analyzed in \cite{Zelditch2}. In \cite{Zelditch3}, \cite{Zelditch4} and \cite{Zelditch5}, the wave invariants associated to bouncing ball orbits are explicitly calculated using Feynmann diagrams to analyze the stationary phase computation from which Balian-Bloch (wave trace) formulas are derived (see also \cite{Zelditch1}). Under mild dynamical conditions and some additional axial symmetry assumptions, these coefficients can be used to determine the Taylor series of a local boundary parametrization. In particular, this allows one to deduce that such domains are spectrally determined amongst a rich class of analytic, symmetric domains.\\
\\
The wave invariants have also proved useful in variational inverse problems, going back to the seminal papers \cite{GK1} and \cite{GK2}, where the authors proved spectral rigidity for closed manifolds with negative sectional curvature. This was recently generalized to Anosov surfaces in \cite{UPS}. In the setting of bounded domains, it is proved in \cite{HeZe12} that ellipses with Dirichlet/Neumann boundary conditions are infinitessimally spectrally rigid through smooth domains with the same symmetries. These results as well as spectral determination of the Robin function in \cite{GuMe79a} were generalized to Liouville billiard tables of classical type in \cite{PopTop03} and \cite{PopTop12}. They were also extended to the Robin setting in \cite{Vig18}, where both the domain and Robin function were allowed to vary simultaneously. A recent breakthrough was obtained in \cite{HeZe19}, where the authors showed that ellipses of small eccentricity are spectrally determined. Thorough surveys of the inverse spectral problem are contained in \cite{ZelditchSurvey2014}, \cite{ZelditchSurvey2}, \cite{HezariDatchevSurvey} and \cite{Melrosesurvey}.
\\
\\
The present article is inspired by \cite{MaMe82}, where the following theorem is proved:
\begin{theo}[\cite{MaMe82}]\label{Melrose theorem}
If $\Omega$ is a bounded and strictly convex planar region, there exists $N = N(\Omega)$ such that if $j > N$, then the contribution $\widehat{\sigma}_j$ of $E_j$ to $\widehat{\sigma}_D$ is of the form
\begin{align}\label{Melrose 1}
\widehat{\sigma}_j = \frac{1}{2\pi} \int_0^\infty \int_{0}^L e^{i (t - \mu_j(s)) \tau } a(s,\tau) ds d\tau,
\end{align}
where in terms of an arclength coordinate $s$ on $\d \Omega$,
\begin{align*}
\mu_j(s) = \Psi_j(s,s')\big|_{s = s'},
\end{align*}
with $\psi(s,s') = L(\mathfrak{g})$, with $\mathfrak{g}$ the length of a $j$-fold geodesic from $s$ to $s'$ and $a_j$ is periodic in $s$, classical and elliptic of order zero, with principal part of the form $e^{i\pi r_j/4} \alpha_j(s)$, $\alpha_j(s) > 0$.
\end{theo}
\noindent Here, $E_j(t)$ is the cosine kernel associated to the parametrix for the $j$ reflection wave operator constructed in \cite{Ch76}, which is reviewed in Section \ref{Chazarain's Parametrix}. In particular, if $L_j$ is a simple length corresponding to a nondegenerate periodic orbit, then Theorem \ref{Melrose theorem} gives an asymptotic expansion for the localized wave trace. In fact, if $\Omega$ satisfies the noncoincidence condition \eqref{NCC}, i.e. $|\d \Omega|$ is not a limit point from below of the lengths of orbits of rotation number $m/n$ for $m \geq 2$, then the trace in Theorem \ref{Melrose theorem} is a spectral invariant. The purpose of this article is to explicitly calculate the principal symbol $a(s,\tau)$ of \eqref{Melrose 1} in terms of geometric data associated to the billiard map, both in the case of simple lengths corresponding to nondegenerate periodic orbits and also for one parameter families of degenerate periodic orbits,  all having the same length associated to a caustic of rational rotation number. It also corrects several errors in literature (including Theorem \ref{Melrose theorem} above) on the order of $a \in S_{\text{cl}}^{1/2}(\d\Omega)$.

\section{Billiards}\label{Billiards}

\noindent Before obtaining a singularity expansion for the wave trace, we first review the relevant background needed on billiards. This will be useful in our discussion of Chazarain's parametrix in Section \ref{Chazarain's Parametrix}. In this section, we denote by $\Omega$ a bounded strictly convex region in $\mathbb{R}^2$ with smooth boundary. This means that the curvature of $\d \Omega$ is a strictly positive function. The billiard map is defined on the coball bundle of the boundary $B^* \d \Omega = \{(q, \zeta) \in T^*\d \Omega : |\zeta| < 1 \}$, which can be identified with the inward part of the circle (cosphere) bundle $S_{\d \Omega}^* \mathbb{R}{^2}$, via the natural orthogonal projection map. We can also identify $B^* \partial \Omega$ with ${\mathbb{R}}\slash {\ell \mathbb{Z}} \times (0, \pi)$, where $\ell=|\partial \Omega|$ is the length of the boundary. Define
\begin{align*}
t_{\pm}^1(y,\eta)& = \inf\{t > 0 : g^{\pm t}(y,\eta) \in  \d \Omega \},\\
t_{\pm}^{-1}(y,\eta) &= \sup \{t < 0 : g^{\pm t}(y,\eta) \in  \d \Omega \},
\end{align*}
where $\pi_1$ is projection onto the first factor and $g^{\pm t}$ is the forwards $(+)$ or backwards $(-)$ geodesic flow on $\mathbb{R}^2$, corresponding to the Hamiltonian $H_\pm = \pm|\eta|$ (see Section \ref{Chazarain's Parametrix}).
If $(y,\zeta) \in B^* \d \Omega$ is mapped to the inward pointing covector $(y,\eta) \in T_{\d \Omega}^* \mathbb{R}^2$ under the inverse projection map, then we define
$$
\beta^{\pm1}(y,\eta) = \widehat{g^{t_\pm^1}(y,\eta)},
$$
where a point $\widehat{(x,\xi)}$ is the reflection of $\xi$ through the cotangent line $T_x^*\d \Omega$. In otherwords, $\widehat{(x,\xi)}$ has the same footpoint and (co)tangential projection as $(x,\xi)$, but reflected conormal component, so that it is again in the inward facing portion of the circle bundle. We call $\beta := \beta^{+1}$ the billiard map. It is well known that $\beta$ preserves the natural symplectic form induced on $B^*(\d \Omega)$ and is differentiable there, extending continuously up to the boundary. The maps $\beta^{\pm n}$ are defined via iteration and it is clear that $\beta^{-n} = (\beta^{n})^{-1}$ for each $n \in \mathbb{Z}$. Associated to the billiard map is the billiard flow, or broken bicharacteristic flow, which we denote by $\Phi^t$.\\
\\
Geometrically, a billiard orbit corresponds to a union of line segments which are called links. A smooth closed curve $\mathcal C$ lying in $\Omega$ is called a {caustic} if any link drawn tangent to $\mathcal C$ remains tangent to $\mathcal C$ after an elastic reflection at the boundary of $\Omega$. By elastic reflection, we mean that the angle of incidence equals the angle of reflection at an impact point on the boundary. We can map $\mathcal C$ onto the total phase space $B^* \partial \Omega$ to obtain a smooth closed curve which is invariant under $\beta$. If the dynamics are integrable, these invariant curves are precisely the Lagrangian tori which folliate phase space. A point $P$ in $B^*\partial \Omega$ is called $q$-periodic ($q \geq 2$) if $\beta^q(P)=P$. We define the rotation number of a $q$-periodic orbit $P$ by $\omega(P)= \frac{p}{q}$, where $p$ is the winding number of the orbit generated by $P$, which we now define. We may consider the modified billiard map $\widetilde{\beta} = \Pi^* \beta$, where $\Pi$ is the natural mapping from $\mathbb{R} / \ell \mathbb{Z} \times [0,\pi]$ to the closure of the coball bundle $\overline{B^*\d \Omega}$. Pulling back by $\Pi$ clearly preserves the notion of periodicity. There exists a unique lift $\widehat{\beta}$ of the map $\widetilde{\beta}$ to the closure of the universal cover $\mathbb{R} \times [0,\pi]$ which is continuous, $\ell$ periodic in first variable and satisfies $\widehat{\beta}(x,0) = (x,0)$. Given this normalization, for any point $(x,\theta) \in \mathbb{R}/\ell \mathbb{Z} \times [0,\pi]$ in a $q$ periodic orbit of $\widetilde{\beta}$, we see that $\widehat{\beta}^q(x,\theta) = (x + p \ell, \theta)$ for some $p \in \mathbb{Z}$. We define this $p$ to be the winding number of the orbit generated by $\Pi (x,\theta) \in \overline{ B^* \d \Omega}$. We see that even if a point $\Pi (x,\theta)$ generates an orbit which is not periodic in the full phase space but is such that $\pi_1 (\widetilde{\beta}^q(x,\theta)) = x$ for some $q \in \mathbb{Z}$, we can still define a winding number in this case. Such orbits are called loops or geodesic loops. For a given periodic orbit, the winding number is independent of which point in the orbit is chosen, so we sometimes write $\omega(\gamma) = \omega(P)$ for any $P \in \gamma = \{P, \beta(P), \cdots, \beta^{q-1}(P)\}$. For deeper results and a more thorough introduction to dynamical billiards, we refer the reader to \cite{Tabachnikov}, \cite{Katok}, \cite{Popov1994} and \cite{PopovTopalov}.
\\
\\
What will be crucial for us in later sections is a description of all orbits making a fixed number of reflections which connect interior points near the diagonal of the boundary in approximately one rotation. These orbits will allow us to cook up phase functions in Section \ref{Chazarain's Parametrix} which parametrize the canonical relation of the wave propagator.
\begin{theo}[8 Orbit Theorem]{\label{8 orbit lemma}}
There exist $C_0 > 0$ and $j_0 = j_0(\Omega)$ sufficiently large such that for $j \geq j_0$ and any two points $x,y \in \text{int}(\Omega)$ which are $C_0/j^4$ close to the diagonal of the boundary, there exist precisely four distinct, broken geodesics of $j$ reflections making approximately one counterclockwise rotation, emanating from $x$ and terminating at $y$. Similarly, there exist four such orbits in the clockwise direction.
\end{theo}

\noindent The proof of Theorem \ref{8 orbit lemma} is based on several lemmas in Section \ref{Proof of 8 orbit lemma} below and is inspired by the author's previous work in \cite{Vig18}, where a similar construction is adapted to elliptical billiard tables. As in that paper, the proof actually provides more information. We now explain what is meant by approximately one rotation. Let $\xi \in S_x^* \Omega$ be one of the $4$ covectors corresponding to the initial condition of a counterclockwise orbit described in Theorem \ref{8 orbit lemma}. Denote by $\widehat{x} = \pi_1 g^{t_1^+}(x,\xi)$ the first point of reflection at the boundary ($\pi_1$ is projection onto the first factor) and by $\widehat{y}$ the $(j+1)$st point of reflection at the boundary after the orbit reaches $y$. If $x, y$ are $O(j^{-1})$ close to the diagonal of the boundary, then $|\widehat{x} - \widehat{y}| = O(j^{-1})$ (see Section \ref{Proof of 8 orbit lemma}). Also let $\omega$ be the angle of reflection made by the orbit at $\widehat{x}$ and note that $\widehat{x}, \widehat{y}$ and $\omega$ all depend implicitly on $\xi$. By approximately one counterclockwise rotation, we mean that for each of the initial covectors $\xi \in S_x^*\Omega$ of the $4$ counterclockwise orbits provided by Theorem \ref{8 orbit lemma}, we have
\begin{equation*}
|\pi_1 \widehat{\beta}^j(\widehat{x}, \omega) - \widehat{y} - \ell| \leq \ell/100.
\end{equation*}
\noindent Here, $\ell = |\d \Omega|$ and $\widehat{\beta}$ is the lift of the billiard map to the closure of the universal cover $\mathbb{R} \times [0,\pi]$ as described in Section \ref{Proof of 8 orbit lemma}. The choice of $\ell/100$ is somewhat arbitrary, but having $\ell$ in the numerator allows for scale invariance and finding the optimal constant is irrelevant for our purposes. The notion of approximately one clockwise rotation is defined analogously.

\begin{def1}\label{TTNN}
Of the four counterclockwise orbits emanating from $x$, two of them become tangent to a level curve of the distance function $d(z) = \text{dist}(z, \d \Omega)$ before making a reflection at the boundary. We denote these orbits by $T$ orbits (for tangency) and call their first links $T$ links. The other two orbits make a reflection at the boundary before becoming tangent to a level curve of $d$ and we call these $N$ orbits (for nontangency); their first link is called an $N$ link. Within either $T$ or $N$ category for the first link, the final link of one of the orbits reaches $y$ before becoming tangent to a level curve (an $N$ link) and the other has a point of tangency before reaching $y$ (a $T$ link). In this way, we obtain four types of counterlockwise orbits from $x$ to $y$, which we denote by $TT$, $TN$, $NT$, and $NN$. The same classification also applies to the four clockwise orbits.
\end{def1}
\begin{figure}
\includegraphics[scale=0.25]{TTorbit1.png}
\includegraphics[scale=0.25]{TNorbit1.png}
\includegraphics[scale=0.25]{NTorbit1.png}
\includegraphics[scale=0.25]{NNorbit2.png}
\caption{Counterclockwise orbit configurations TT, TN, NT, and NN corresponding to $j = 4$. The green and pink curves are the distance curves on which $x$ and $y$ respectively lie. The red curve is a caustic to which the billiard orbit is tangent in the completely integrable setting.}
\label{configurations}
\end{figure}
\noindent See Figure \ref{configurations} for an example on the ellipse with $j=4$. These configurations will be important in determining which limiting orbits give rise to geodesic loops of precisely $j$ reflections as $(x,y) \to \Delta \d \Omega$, where $\Delta: \d \Omega \to \d \Omega \times \d \Omega$ is the diagonal embedding $x \mapsto (x,x)$.
\begin{def1}\label{psi jk}
For $1 \leq k \leq 8$, we set $\Psi_j^k(x,y)$ to be a branch of the length functional corresponding to one of the orbits in Theorem \ref{8 orbit lemma}. It depends only on $x,y, j$ and $k$. We use the convention that for a fixed number of reflections $j$, the indices $1 \leq k \leq 4$ correspond to the counterclockwise orbits $TT, TN, NT, NN$ in that order and the indices $5 \leq k \leq 8$ correspond to their clockwise counterparts in the same order.
\end{def1}
\noindent The author learned of a similar function in \cite{MaMe82} (page 492), where its restriction to the boundary is defined. In such a case, i.e. if $x, y \in \d \Omega$, it is stated in \cite{MaMe82} and proved in \cite{Popov1994} that only a single counterclockwise orbit of $j$ reflections exists between the boundary points if they are sufficiently close and $j$ is sufficiently large.  The proof in Section \ref{Proof of 8 orbit lemma} below also shows that as $x$ and $y$ approach the diagonal of the boundary from the interior of $\Omega$, the corresponding orbits coalesce and converge to the orbits described in \cite{MaMe82}. This in fact proves the claims made in \cite{MaMe82} and provides a simpler alternative to the methods employed in \cite{Popov1994}. The limiting orbits may have a different number of reflections, which is addressed in the proof of Lemma \ref{Sj lemma}.

\section{Proof of Theorem \ref{8 orbit lemma}}\label{Proof of 8 orbit lemma}

\noindent We will construct conic subbundles $C_\pm^*(\Omega; j)$ of $T^*\Omega$ over a tubular neighborhood $\widetilde{U}_j$ of the boundary with the following property: for each $x_0 \in \widetilde{U}_j$, there exists a conic subset of $T_{x_0}^*\Omega \backslash 0$ such that all orbits emanating from $x_0$ which make $j$ reflections and approximately one rotation in the sense of Section \ref{Billiards} have initial covector in $C_{x_0}^*(\Omega;j)$, the fiber of $\cup_\pm C_\pm^*(\Omega;j)$ at $x_0$. At each $x_0 \in \widetilde{U}_j$, there also exist distinguished covectors $\xi_0^\pm \in T_{x_0}^*(\Omega) \backslash 0$ which are tangent to distance curves folliating a neighborhood of the boundary. If $\widetilde{U}_j$ is chosen sufficiently small, the $j$ reflection orbits emanating from $\xi_0^\pm$ will make less than a quarter rotation. By rotating $\xi$ away from $\xi_0^\pm$ in either direction within the fiber of $C_\pm^*(\Omega;j)$ at $x_0$, we will show that the angle of reflection at the first impact point on the boundary increases monotonically. From the boundary, we can then take advantage of the twist property of the billiard map, and show that monotonicity of the incident angles causes orbits of a large number of reflections to wind around $\Omega$. Two of the orbits from $x$ to $y$ will be obtained by perturbing $\xi_0^+$ in the counterclockwise direction and the other two will be obtained by perturbing $\xi_0^+$ in the clockwise direction. The four clockwise orbits will then be constructed in a similar manner, rotating $\xi_0^-$ in both the clockwise and counterclockwise fiber directions. The arguments in this section will also provide the additional topological structure of the resulting orbits referenced in Definition \ref{TTNN}, which can be seen in Figure \ref{configurations}.

\subsection{Lazutkin coordinates} Recall that billiard map is defined on the coball bundle $B^*\d \Omega$, which can be identified with the collection of inward facing covectors in the circle bundle $S_{\d \Omega}^{*} \mathbb{R}^2$. Letting $s$ denote the arclength parameter on $\d \Omega$ and $\vartheta$ the angle an inward facing covector makes with the positively oriented boundary, we defined the modified billiard map $\widetilde{\beta}$ in arclength coordinates $(s, \vartheta) \in \mathbb{R}/\ell \mathbb{Z} \times (0,\pi)$. In this coordinate system, the modified billiard map is given by $\widetilde{\beta}(s_1, \vartheta_1) = (s_2, \vartheta_2)$, where
\begin{align}
\begin{cases}
s_2 = s_1 + a_1(s_1) \vartheta_1 + a_2(s_1) \vartheta_1^2 + a_3(s_1) \vartheta_1^3 + F(s_1, \vartheta_1)\vartheta_1^4\\
\vartheta_2 = \vartheta_1 + b_2(s_1) \vartheta_1^2 + b_3(s_1) \vartheta_1^3 +  G(s_1, \vartheta_1)\vartheta_1^4,
\end{cases}
\end{align}
and $a_i, b_i, F$ and $G$ are smooth functions. This is a Taylor expansion in the angular variable near $\vartheta_1 = 0$. There are explicit formulas for the coefficients $a_i, b_i$. In particular, $a_1(s_1) = 2 \rho(s_1)$ and $b_2 = -2/3 \rho'(s_1)$, where $\rho(s_1)$ is the radius of curvature at $s_1$ (see \cite{Lazutkin}). In \cite{Lazutkin}, the change of coordinates
\begin{align}\label{Lazutkin coord}
x = \frac{\int_0^{s} \rho^{2/3}(t)dt}{\int_0^{\ell} \rho^{2/3}(t)dt}, \qquad
\alpha = \frac{4 \rho(s) \sin \vartheta/2}{\int_0^{\ell} \rho^{2/3}(t)dt}
\end{align}
was introduced near the circle corresponding $\vartheta = 0$. We call the coordinate system $(x, \alpha)$ in \eqref{Lazutkin coord} Lazutkin coordinates. The advantage of this change of variables is that in these coordinates, the billiard map becomes a small perturbation of the translation map
\begin{align}
\begin{pmatrix}
x_1\\
\alpha_1
\end{pmatrix} \mapsto
\begin{pmatrix}
1 & 1\\
0 & 1
\end{pmatrix}
\begin{pmatrix}
x_1\\
\alpha_1
\end{pmatrix}
= \begin{pmatrix}
x_1 + \alpha_1\\
\alpha_1
\end{pmatrix},
\end{align}
which is completely integrable. The billiard map is given in Lazutkin coordinates by $(x_1, \alpha_1) \mapsto (x_2, \alpha_2)$, where
\begin{align}
\begin{cases}
x_2 = x_1 + \alpha_1 + \alpha_1^3 f(x_1, \alpha_1)\\
\alpha_2 = \alpha_1 + \alpha_1^4 g(x_1,\alpha_1),
\end{cases}
\end{align}
for some smooth functions $f,g$. The original implementation of these coordinates was in \cite{Lazutkin}, where it was shown as a consequence of the KAM theorem (see \cite{KAM1}, \cite{KAM4}, \cite{KAM2}, \cite{KAM3}) that there exist an abundance of invariant tori for the billiard map near $\alpha = 0$. However, we will only use these coordinates to prove Theorem \ref{8 orbit lemma}. Without loss of generality, we will often interchange the use of arclength and Lazutkin coordinates in the domain of the billiard map $\beta$.

\subsection{Angles of reflection}
Before estimating the billiard map and its iterates from the boundary, we need several preliminary estimates on the angles of reflection. In this section, we denote $\widehat{\beta}^k(x_1, \alpha_1) = (x_k, \alpha_k)$.
\begin{lemm}\label{angle bound 1}
For each $c_1, c_2 > 0$, there exists $j_1 = j_1(\Omega, c_1, c_2) \in \mathbb{N}$ sufficiently large such that for $j \geq j_1$ and $c_1/j \leq \alpha_1 \leq c_2/j$, it follows that
$$
\frac{c_1}{(1 + c_1/j)^k j} \leq \alpha_k = \pi_2 \beta^k(x_1,\alpha_1) \leq \frac{c_2 (1+ \frac{c_2}{j})^k}{j} \qquad (1 \leq k \leq j),
$$
where $\pi_2$ is the projection onto the angular variable. In particular, $c_1 e^{-c_1}/j \leq \pi_2 \beta^k(x_1,\alpha_1) \leq c_2 e^{c_2}/j$.
\end{lemm}
\begin{proof}
We use Lazutkin coordinates and proceed by induction. Supose the claim is true for $1, \cdots, j - 1$. Then
\begin{align}
\begin{split}
\pi_2 \beta^j (x_1,\alpha_1) &= \pi_2 \beta \circ \beta^{j-1} (x_1, \alpha_1) = \alpha_{j-1} + \alpha_{j-1}^4 g(x_1,\alpha_1)\\
&\leq c_2(1 + c_2/j)^{j-1}/j + (c_2(1+ c_2/j)^{j-1}/j)^4 g(x_1,\alpha_1)\\
&= \frac{c_2(1 + c_2/j)^{j-1}}{j} \left(1 + \frac{c_2^3(1 + c_2/j)^{3j-3} g(x_1,\alpha_1)}{j^3} \right).
\end{split}
\end{align}
Choosing $j_1$ sufficiently large so that
\begin{align}
\frac{c_2^3(1+c_2/j_0)^{3j_0 - 3}g(x_1,\alpha_1)}{j_0^3} \leq \frac{c_2}{j_0},
\end{align}
the upper bound follows. Similarly,
\begin{align}
\begin{split}
\pi_2 \beta^j (x_1,\alpha_1) &= \pi_2 \beta \circ \beta^{j-1} (x_1, \alpha_1) = \alpha_{j-1} + \alpha_{j-1}^4 g(x_1,\alpha_1)\\
&\geq \frac{c_1}{(1 + c_1/j)^{j-1} j} + \left(\frac{c_1}{(1+ c_1/j)^{j-1}j}\right)^4 g(x_1,\alpha_1)\\
&= \frac{c_1}{(1 + c_1/j)^{j-1}j} \left(1 + \frac{c_1^3}{(1 + c_1/j)^{3j-3} j^3} g(x_1,\alpha_1) \right),
\end{split}
\end{align} which can be made greater than
\begin{align*}
\frac{c_1}{(1 + c_1/j)^k j},
\end{align*}
under the condition that
\begin{align*}
\frac{c_1}{j_1} + \frac{c_1^3 g(x_1,\alpha_1)(1 + c_1/j_1)}{(1+ c_1/j_1)^{3j_1 - 3} j_1^3} \geq 0.
\end{align*}
This certainly holds for $j_1$ large enough, which concludes the proof.
\end{proof}

\begin{lemm}\label{L3}
There exist $C_1 , C_2 > 0$ and $j_2 = j_2(\Omega) \in \mathbb{N}$ with the following property: for all  $j \geq j_2$ and $x, y$ which are $O(1/j)$ close to the diagonal of the boundary, any orbit $\gamma$ of $j$ reflections which emanates from $x$ and terminates at $y$ in approximately one rotation (in the sense of Section \ref{Billiards}) has angles of reflection $\alpha_k$ in the range
$$
\frac{C_1}{j} \leq \alpha_k \leq \frac{C_2}{j} \qquad (1 \leq k \leq j).
$$
\end{lemm}
\begin{proof}
Let $\widehat{x}, \widehat{y}$ denote the first and $(j+1)$st points of reflection on the boundary. Recall that by approximately one rotation, we mean that
\begin{align*}
|\widehat{\beta}^j(\widehat{x}, \alpha_1) - \widehat{y} - \ell| \leq \ell /100,
\end{align*}
where $\widehat{\beta}$ is the lift to the universal cover and $\ell = |\d \Omega|$. Each arc separated by moments of reflection $x_p, x_{p+1}$ on the boundary has length $x_{p+1} - x_p$. As $\gamma$ makes approximately one rotation, we see that there must exist one arc of length $x_{p+1} - x_p \geq (99/100)\ell/j$. Similarly, there must exist an arc of length $x_{m+1} - x_m \leq (101/100) \ell/j$. By Proposition 14.1 of \cite{Lazutkin5}, for each $1 \leq k \leq j$,
\begin{align}\label{curvature bound}
2 \alpha_k \rho_\text{min} \leq x_{k+1} - x_k \leq 2 \alpha_k \rho_\text{max},
\end{align}
where $\rho_\text{min}$ and $\rho_\text{max}$ are the minimum and maximum radii of curvature respectively for $\d \Omega$. In particular,
\begin{align}
\frac{99 \ell }{200 \rho_\text{max} j} \leq \alpha_p, \qquad \alpha_m \leq \frac{101 \ell}{200 \rho_\text{min} j}.
\end{align}
We now apply Lemma \ref{angle bound 1} to $\beta^\pm$ at the points $x_p$ and $x_m$ to conclude that for all $1 \leq k \leq j$,
\begin{align}
\frac{\exp\left( - \frac{99 \ell }{200 \rho_\text{max}}\right)}{j} \leq \alpha_k \leq \frac{\exp\left( \frac{101 \ell}{200  \rho_\text{min}} \right)}{j}.
\end{align}
Here we have used Lemma \ref{angle bound 1} applied to both the billiard map and its inverse, with the constants $c_1 = 99 \ell / 200 \rho_\text{max}$ and $c_2 = 101 \ell /200 \rho_\text{min}$ which depend only on $\Omega$. Hence, $j_2$ may be chosen uniformly for all orbits making approximately one rotation regardless of their initial positions.
\end{proof}

\noindent We now investigate the angle of reflection at the first point of impact on the boundary. If $d(z) = \text{dist}(z, \d \Omega)$, then the level curves of $d$ folliate a neighborhood of the boundary. For $x_0$ in the tubular neighborhood $\widetilde{U}_j$ (whose diameter remains to be specified), we define $\xi_0^\pm \in \pm T_{x_0}^* d^{-1}(d(x_0)) \cap S_{x_0}^*\Omega$ so that $\xi_0^+$ points in the counterclockwise direction and $\xi_0^-$ points in the clockwise direction.
\begin{figure}
\begin{center}
\begin{tikzpicture} [domain=0:2, scale = 1.1]
\draw[->, thick] (-5,0)--(5,0);
\draw[->, thick] (0,-2)--(0,4);
\draw[->, thick] (0,-2)--(0,4);
\draw[black, line width = 0.25mm]   plot[smooth,domain=-4:4] (\x, 0.2*\x*\x);
\coordinate (A) at (0, 2.5);
\coordinate (B) at (2,0.8);
\coordinate (C) at (0.25,-0.6);
\coordinate (D) at (4, 2.4);
\coordinate (E) at (-3.53553390593,2.5);
\coordinate (F) at (3.53553390593,2.5);
\coordinate (a-) at (-3.53553390593, 0);
\coordinate (a+) at (3.53553390593, 0);
\draw[black, thick] (A)--(B);
\draw[black, thick] (C)--(D);
\draw[black, thick, dashed] (E)--(F);
\draw[black, thick, dashed] (a-)--(E);
\draw[black, thick, dashed] (a+)--(F);
\draw[fill=black] (a-) circle (2pt) node[below left]{$a_-$};
\draw[fill=black] (a+) circle (2pt) node[below right]{$a_+$};
\draw[fill=black] (A) circle (2pt) node[above left]{$b$};
\draw[fill=black] (B) circle (2pt) node[below right]{$(x,f(x))$};
\path
  pic["$\alpha_1$",draw=black,<->,angle eccentricity=1.2,angle radius=0.75cm] {angle=A--B--C};
\path
  pic["$\theta$",draw=black,<->,angle eccentricity=1.2,angle radius=0.75cm] {angle=B--A--F};
\end{tikzpicture}
\end{center}
\caption{Angular derivative of the billiard flow.}
\label{picture of billiard flow}
\end{figure}
Let $\xi \in S_{x_0}^* \Omega$ be identified with the corresponding point $\theta \in S^1 = \mathbb{R}/\mathbb{Z}$, which denotes the clockwise angular parametrization of the fiber $S_{x_0}^*\Omega$, normalized so that $\xi_0^+$ and $\xi_0^-$ correspond to $0$ and $\pi$ respectively. Denote by $\alpha_1 \in (0,\pi)$ the angle made between the billiard orbit emanating from $(x_0, \xi)$ and the positively oriented cotangent line at the first point of reflection at the boundary. It is clear that $|\d \xi / \d \theta| = 1$.
\begin{lemm}\label{L1}
There exist $c = c(\Omega) > 0$ and $j_3 = j_3(\Omega) \in \mathbb{N}$ such that if $j \geq j_3$ and $d(x_0) = O(1/j^2)$, then
$$
\frac{\d \alpha_1}{\d {\theta}} \geq c
$$
for all initial conditions $(x_0, \xi(\theta))$ which generate a counterclockwise orbit of $j \geq j_0$ reflections, making approximately one rotation.
\end{lemm}

\begin{proof}
Let $x_0 \in \text{int} (\Omega)$ be near the boundary and consider the point ${x} \in \d \Omega$ which minimizes the Euclidean distance $|x_0- x|$ among all boundary points. We then apply an affine change of coordinates $y \mapsto A(y - x)$ with $A \in SO(2)$ so that $x$ is mapped to the origin and the positively oriented unit tangent vector at ${x}$ is mapped to $(1,0)$. The vector ${x} - x_0$ is perpendicular to both $T_x^*\d \Omega$ and $T_{x_0}^* d^{-1}(|x - x_0|)$ (see Figure \ref{picture of billiard flow}). Therefore, the point $A(x_0 - x)$ lies on the positive vertical axis and we denote it by $(0,b)$. In these coordinates, the boundary $\d \Omega$ is locally parametrized by the graph of a smooth convex function $f$. Since $f'(0) = 0$, we see that $f(t) = \kappa t^2 + R(t)$, with $\kappa$ denoting the curvature at the point corresponding to $(0,0)$ and $R(t) = O(|t|^3)$. Dilating by a factor of $\kappa^{-1}$, we see that in rescaled coordinates, the function becomes $f(t) =  t^2 + R(t)$, with $R(t)$ again of order $O(|t|^3)$.
\\
\\
In graph coordinates, $\alpha_1$ becomes the angle between $(t, f(t)) - (0,b)$ and the negatively oriented cotangent line $T^*_{(t,f(t))} \d \Omega$ as illustrated in Figure \ref{picture of billiard flow}. This shows that
\begin{align}\label{cos angle}
\cos \alpha_1 = \frac{(t,f(t) - b) \cdot (1, f'(t))}{\sqrt{t^2 + (f(t) - b)^2} \sqrt{1+ (f'(t))^2}} \geq 0,
\end{align}
since $0 \leq \alpha_1 \leq \pi/2$. Recall that $\xi_0^\pm$ are the unit covectors in the positive (+), resp. negative (-), cotangent line $T^*_{x_0} d^{-1}(|x- x_0|)$ corresponding to $\theta = 0$ and $\theta = \pi$ respectively. For now, we only consider counterclockwise orbits obtained by perturbing the orbit emanating from $(x_0,\xi_0^+)$, as they correspond to $-\pi/2 \leq \theta \leq \pi/2$ in the statement of the lemma. Denote the points $a_- = \min\{f^{-1}(b)\}$ and $a_+ = \max \{f^{-1}(b)\}$. We first show that $d \alpha_1/dt < 0$ on $[0,a_+]$ so that the angle of incidence at the first impact point on the boundary is increasing as the initial covector of the trajectory winds clockwise in $S_{x_0}^*\Omega$ away from $\xi_0^+$, i.e. $d\alpha/dt < 0$ on $[0,a_+]$.
\\
\\
Noting that $t$ and $\theta$ are negatively correlated, it suffices to show that the logarithmic $t$ derivative of $\cos \alpha_1(t)$ is positive:
\begin{align}\label{log cos derivative}
\frac{d}{dt} \log \cos \alpha_1(t) = \frac{(1 + f'' (f - b) + (f')^2)}{t + f' (f-b)} - \frac{t + f'(f-b)}{t^2 + (f-b)^2} - \frac{f'f''}{1 + (f')^2}.
\end{align}
Multiplying \eqref{log cos derivative} through by a common denominator, which is positive for $t, b$ small, we obtain
\begin{align}\label{logarithmic derivative times common den}
\begin{split}
g(t,b) :=& (t + f'(f - b))(t^2 + (f - b)^2)(1 + (f')^2) \frac{d}{dt} \log \cos \alpha(x)\\
=& (1 + f'' (f - b) + (f')^2)(t^2 + (f - b)^2)(1 + (f')^2)\\
&- (t + f'(f - b))^2(1 + (f')^2)  - f' f''(t +f'(f-b))(t^2 + (f-b)^2).
\end{split}
\end{align}
To show that \eqref{logarithmic derivative times common den} is positive, we plug in our second order Taylor approximation for $f$ and expand $g$ to fourth order in a parabolic neighborhood of $t= 0, b = 0$:
\begin{gather}\label{fourth order expansion}
g(t,b) = 3t^4 + 10t^2 b^2 -2b^3 + b^2 + O((t^2 + b^2)^{5/2}).
\end{gather}
As $f(t) = t^2 + O(t^3)$, we see that $a_\pm = \pm b^{1/2} + O(b^{3/2})$, so we choose $b \leq 1/j^2$ sufficiently small in order to make \eqref{fourth order expansion} positive in a parabolic neighborhood of origin. Examining the initial set up, choosing $b$ small amounts to choosing $x_0$ close to the boundary, all of which can be done uniformly in $x_0$ and the curvature $\kappa$. Choosing a uniform $b$ gives us a tubular neighborhood of the boundary $\widetilde{U}_j = \{x: 0< d(x,\d \Omega) < j^{-2} \}$, as mentioned in the begining of the section. We will later shrink $\widetilde{U}_j$ by a factor of $j^{-2}$, following Lemma \ref{clockwise int points}.
\\
\\
While the lower bound in \eqref{fourth order expansion} is of order $b^2 = O(j^{-4})$, we in fact need uniform bounds for $\d \alpha_1/ \d \theta$ which we now provide. The prefactor in \eqref{logarithmic derivative times common den} is nonvanishing for $t \neq 0, t \neq a_+$ and in particular, for all $t$ corresponding to an orbit which makes approximately one rotation in $j$ reflections. We first need to compare $\d/\d t$ and $\d / \d \theta$ in terms of $b$, which depends on $j$. In graph coordinates, $\theta$ becomes the angle of the first link from the horozontial axis (see Figure \ref{picture of billiard flow}). Hence, $\theta = \arctan((b- f(t))/t)$ and
\begin{align*}
\frac{\d \theta}{\d t} = \frac{1}{1 + (b - f(t))^2/t^2 } \left( - \frac{ f'(t)}{t} - \frac{b - f(t)}{t^2}\right),
\end{align*}
or equivalently,
\begin{align}\label{d x d xi}
\frac{\d t}{\d \theta} = - \left(\frac{t^2 + (b - f(t))^2}{t f'(t) + b - f(t)}\right),
\end{align}
where $t$ is defined implicitly as a function of $\theta$. Plugging \eqref{d x d xi} into  \eqref{log cos derivative} and using that \eqref{fourth order expansion} is bounded below by $b^2/2$ for $b$ sufficiently small, we see that
\begin{align}\label{dalphadxi}
\begin{split}
\frac{\d \alpha}{\d \theta} &\geq \frac{ b^2 (t^2 + (b - f(t))^2) }{2 (tf'(t) + b - f(t)) (t + f'(t)(f(t) - b)) (t^2 + (f(t) - b)^2 )(1 + (f'(t))^2)  \tan \alpha_1  }\\
&=  \frac{ b^2 }{2 (tf'(t) + b - f(t) ) (t + f'(t)(f(t) - b)) (1 + (f'(t))^2)  \tan \alpha_1  }.
\end{split}
\end{align}
The denominator in \eqref{dalphadxi} is nonvanishing and can be estimated above on $(0, a_+]$ by
\begin{align}
c' t b \tan \alpha_1 + O(b^{5/2}).
\end{align}
for some $c' > 0$, independent of $j$. Observe that if $t \leq b$ for example, then the denominator is bounded above by $c' b^2$ and $\d \alpha_1/ \d \theta \geq c''$, for some $c'' > 0$ which is also independent of $j$.
\\
\\
We now find a similar upper bound on $t(\theta)$ in terms of $b$, corresponding to the set of covectors $\xi(\theta)$ which produce orbits of $j$ reflections and approximately one rotation. Let $b = r b_0$, where $b_0 = j^{-2}$ for $r \in (0,1]$. Recall from Lemma \ref{L3} that there exist constants $C_1, C_2 > 0$ such that for all orbits making approximately one rotation and $j$ reflections, we have $C_1/j \leq \alpha_k \leq C_2/j$ for each $1 \leq k \leq j$. Hence,
\begin{align}\label{b_0 bound}
C_1 b_0^{1/2} \leq \alpha_1 \leq C_2 b_0^{1/2}.
\end{align}
Observe that equation \eqref{cos angle} gives
\begin{align}\label{log cos alpha}
\begin{split}
\log \cos \alpha_1(t) &= \log(t + f'(t)(f(t) - b)) - 1/2 \log (t^2 + (f(t)-b)^2)\\
&- 1/2 \log(1 +(f'(t))^2).
\end{split}
\end{align}
Denote the terms in \eqref{log cos alpha} above by $A, B$ and $C$. Then,
\begin{align*}
\begin{split}
A &= \log( t + (2t + O(t^2) ) (t^2 + O(t^3) - b ))\\
&= \log t + \log(1 + 2t^2 - 2b + O(t^3) + O(b))\\
&= \log t + 2t^2 + O(t^3) + O(b).
\end{split}
\end{align*}
Similarly,
\begin{align*}
\begin{split}
B &= -1/2 \log(t^2 + t^4 - 2t^2b + b^2 + O(t^5) + O(t^3 b) )\\
&= -\log t - 1/2 (t^2 - 2b + b^2/t^2 + O(t^3) + O(tb))
\end{split}
\end{align*}
and
\begin{align*}
C = - \log(1 + (2t + O(t^2)) ) = -4t^2 + O(t^3).
\end{align*}
As $t = O(b^{1/2})$, we have
\begin{align*}
A + B + C = \frac{-b^2}{2t^2} + O(b).
\end{align*}
On the other hand,
\begin{align*}
\begin{split}
A + B + C &= \log \cos \alpha_1 = \log (1 - a_1^2 / 2 + O(\alpha_1^4)).\\
&= - \alpha_1^2/2 +  O(b_0^{4/N}),
\end{split}
\end{align*}
which implies that
\begin{align*}
b^2/t^2 = r^2 b_0^2 / t^2 \sim \alpha_1^2 / 2.
\end{align*}
In particular, \eqref{b_0 bound} implies
\begin{align}\label{x sim}
t \sim \frac{\sqrt{2}r b_0}{\alpha_1} \in \left[\frac{\sqrt{2} rb_0^{1/2}}{C_2} , \frac{\sqrt{2} rb_0^{1/2}}{C_1} \right].
\end{align}
Note also that $\tan \alpha_1 \leq 2 \alpha_1 \leq 2 C_2 b_0^{1/2}$. Combining this with \eqref{x sim}, we see that the denominator in equation \eqref{dalphadxi} is bounded above by
\begin{align}
c' t b \tan \alpha_1 + O(b^{5/2}) \leq c''' b  (r b_0^{1 - 1/2}) b_0^{1/2} = c''' b^2,
\end{align}
for some $c''' > 0$ independent of $j$. This establishes the claim that $\d \alpha/ \d \theta$ is uniformly bounded below on the set of initial covectors corresponding to orbits of $j$ reflections and approximately one rotation.
\end{proof}

\subsection{Monotonicity and the twist property}
The significance of increasing the angle of reflection at the first impact point on the boundary is that after $j-1$ reflections at the boundary, the $jth$ impact point winds around $\d \Omega$ with approximate speed $j$. This is essentially the twist property of the billiard map (see \cite{Tabachnikov} for a general definition of twist maps). The following lemma makes this notion quantitative.
\begin{lemm}\label{boundary monotonicity}
There exists $j_4 = j_4(\Omega) \in \mathbb{N}$ such that if a trajectory makes $j \geq j_4$ reflections at the boundary and approximately one rotation, then the Jacobian of the iterated billiard map satisfies
$$
D \widehat{\beta}^j (x_1 ,\alpha_1) = \begin{pmatrix}
1 & j\\
0 & 1
\end{pmatrix}
+ O(1/j)
$$
in Lazutkin coordinates lifted to the universal cover. In particular, there exist constants $C_3, C_4, C_5, C_6 > 0$ depending only on $\Omega$ such that for any $(x_1, \alpha_1)$, we have
\begin{align*}
C_3 j \leq \frac{\d}{\d \alpha_1} \pi_1 \widehat{\beta}^j(x_1, \alpha_1) \leq C_4 j\\
C_5 \leq \frac{\d}{\d \alpha_1 } \pi_2 \widehat{\beta}^j(x_1, \alpha_1)  \leq C_6,
\end{align*}
where $\pi_1$ and $\pi_2$ denote projections onto the first and second components respectively.
\end{lemm}
\begin{proof}
In Lazutkin coordinates, $\beta(x_k,\alpha_k) = (x_k + \alpha_k + \alpha_k^3 f(x_k,\alpha_k), \alpha_k + \alpha_k^4 g(x_k,\alpha_k))$. Hence,
\begin{align}\label{Dbetaj}
D\widehat{\beta}^j(x_1,\alpha_1) = \prod_{k = 1}^j \begin{pmatrix}
1 + \alpha_k^3 \frac{\d f}{\d x}(x_k, \alpha_k) & 1 + 3\alpha_k^2 f(x_k,\alpha_k) + \alpha_k^3 \frac{\d f}{\d \alpha}(x_k,\alpha_k)\\
\theta_k^4 \frac{\d g}{\d x}(x_k,\alpha_k) & 1 + 4 \alpha_k^3 g(x_k,\alpha_k) + \alpha_k^4 \frac{\d g}{\d \alpha}(x_k, \alpha_k)
\end{pmatrix},
\end{align}
where $(x_k, \alpha_k) = \widehat{\beta}^k(x_1,\alpha_1)$ and $f,g$ are extended to be $\ell$ periodic in $x$ on $\mathbb{R} \times [0,\pi]$. Each of the terms in the product can be decomposed into $(I_2 + B_k)$, where
\begin{gather}\label{Bk}
B_k =  \begin{pmatrix}
\alpha_k^3 \frac{\d f}{\d x}(x_k, \alpha_k) & 1 + 3\alpha_k^2 f(x_k,\alpha_k) + \alpha_k^3 \frac{\d f}{\d \alpha}(x_k,\alpha_k)\\
\alpha_k^4 \frac{\d g}{\d x}(x_k,\alpha_k) & 4 \alpha_k^3 g(x_k,\alpha_k) + \alpha_k^4 \frac{\d g}{\d \alpha}(x_k, \alpha_k)
\end{pmatrix} = \begin{pmatrix}
0 & 1\\
0 & 0
\end{pmatrix} + O(1/j^2).
\end{gather}
Denote the constant matrix on the righthand side above by
\begin{align}
N = \begin{pmatrix}
0 & 1\\
0 & 0
\end{pmatrix}.
\end{align}
We have
\begin{align}
\begin{split}
D \widehat{\beta}^j &= (I_2 + B_k )^j = \sum_{k = 0}^j \binom{j}{k}\begin{pmatrix}
O(j^{-3}) & 1 + O(j^{-2})\\
O(j^{-4}) & O(j^{-3})
\end{pmatrix}^k
\\
&= I_2 + j N + O(1/j)  + \sum_{k = 2}^j \binom{j}{k}(N + O(1/j^2))^k.
\end{split}
\end{align}
Noting that $N$ is nilpotent of order $2$ (i.e. $N^2 = 0$), we can estimate the sum above by
\begin{align}\label{key decomp}
\begin{split}
D\widehat{\beta}^j = \begin{pmatrix}
1 & j\\
0 & 1
\end{pmatrix}
+ O(1/j) + \sum_{k = 2}^j \binom{j}{k} O(j^{-2(k-1)}).
\end{split}
\end{align}
We want to show that the remainder term $\sum_{k = 2}^j \binom{j}{k}  O(j^{-2(k-1)})$ in \eqref{key decomp} is small. Choosing $j$ large so that $O(j^{-2(k-1)}) \leq 2^{-k} / j $,
\begin{align}
\sum_{k \geq 2} O(j^{-2(k-1)}) \leq \frac{1}{j} \sum_{k = 0}^j \binom{j}{k} 2^{-k} = \frac{1}{j}.
\end{align}
Noting that angular derivatives $\d \pi_1 \widehat{\beta}^j / \d \alpha_1$ and $ \d \pi_2 \widehat{\beta}^j / \d \alpha_1$ are the $(1,2)$ and $(2,2)$ components of $D\widehat{\beta}^j$ respectively concludes the proof of the lemma.
\end{proof}

\begin{rema}
As $j$ increases, Lemma \ref{boundary monotonicity} shows that the Poincar\'e map associated to a periodic orbit of rotation number $1/j$ becomes more and more degenerate, corresponding to the accumulation of caustics at the boundary. In fact, there exist higher order Lazutkin coordinates for which the remainder above could be replaced by $O(j^{-N})$ for any $N \in \mathbb{N}$. However, this is not needed in the remainder of the paper.
\end{rema}

\begin{figure}
\begin{center}
\begin{tikzpicture} [domain=0:2, scale = 1.1]
\begin{scope}[rotate =0]
\coordinate (x) at (2,0);
\coordinate (y) at (0,2);
\coordinate (z) at (3.2,-1.2);
\coordinate (a) at (-2.75,2);
\coordinate (b) at (2.75,2);
\coordinate (c) at (5.2,0.133333333);
\coordinate (d) at (1.2, -2.5333333333);
\draw[black, thick] (x)--(z);
\draw[black, thick] (x)--(y);
\draw[black, thick] (a)--(b) node[above]{$T_{x_1}^*\d \Omega$};
\draw[black, thick] (c)--(d) node[below right]{$T_{x_{-1}}^*\d \Omega$};
\draw[black, thick, ->] (x)--(2,0.9) node[above right]{$\xi_0^+$};
\draw[black, thick] (0,0) ellipse (4cm and 2cm);
\draw[black, dashed] (0,0) ellipse (2cm and 1cm);
\draw[black, dashed] (x) circle (0.9cm);
\draw[fill=black] (x) circle (2pt) node[label = {[xshift = 0.7cm, yshift = -0.42cm] $x_0$}]{};
\draw[fill=black] (z) circle (2pt) node[below right]{$x_{-1}$};
\draw[fill=black] (y) circle (2pt) node[above left]{$x_1$};
\fill[fill opacity = 0.8, fill=green!20!white]
    (2,0) -- (2 + 0.9 * cos 110, 0.9*sin 110 ) arc (110:160:0.9cm) -- cycle;
\fill[fill opacity = 0.8, fill=green!20!white]
    (2,0) -- (2 + 0.9 * cos 200, 0.9*sin 200 ) arc (200:250:0.9cm) -- cycle;
\fill[fill opacity = 0.8, fill=blue!20!white]
    (2,0) -- (2 + 0.9 * cos 290, 0.9*sin 290 ) arc (290:340:0.9cm) -- cycle;
\fill[fill opacity = 0.8, fill=blue!20!white]
    (2,0) -- (2 + 0.9 * cos 20, 0.9*sin 20 ) arc (20:70:0.9cm) -- cycle;
\path
  pic["$\alpha_1$",draw=black,<->,angle eccentricity=1.5,angle radius=0.75cm] {angle=x--y--b};
\path
  pic["$\alpha_{-1}$",draw=black,<->,angle eccentricity=1.5,angle radius=0.5cm] {angle=c--z--x};
\end{scope}
\end{tikzpicture}
\end{center}
\caption{A counterclockwise perturbation of $\xi_0^+$. The green cones at $x_0$ are the length one covectors in $C_{x_0,+}^*(\Omega;j)$ and the blue cones at $x_0$ are the length one covectors in $C_{x_0,-}^*(\Omega;j)$.}
\label{picture of counterclockwise perturbation}
\end{figure}
\noindent The point $x_j$ can be obtained in two ways. The first way is by flowing forwards to the boundary point $x_1 := \pi_1 g^{t_1^+}(x_0, \xi)$, iterating the billiard map $j$ times and projecting onto the first component. The second way is by flowing backwards to the boundary point $x_{-1} := g^{t_1^-}(x_0, \xi)$ and iterating the billiard map $j+1$ times. If $\xi$ is rotated in the clockwise direction away from $\xi_0^+$, it is convenient to use the angle $\alpha_{1}$ at $x_{1}$, since it was shown in Lemma \ref{L1} that $\alpha_{1}$ is increasing in $\theta$. If $\xi$ is rotated counterclockwise away from $\xi_0^+$, it is instead advantageous to use the point $(x_{-1}, \alpha_{-1})$, since this again corresponds to a \textit{clockwise} perturbation of $\xi_0^-$. In either regime $\theta < \theta_0^+$ or $\theta > \theta_0^+$, we can then take advantage of the monotonicity of $x_j$ in $\alpha_{\pm 1}$ as $\xi$, or equivalently $\theta$, varies. The next lemma shows that the point $x_j$ is monotonically increasing for all orbits making approximately one rotation as $\xi(\theta)$ winds in either direction.

\begin{lemm}\label{winding in xi}
There exist $C_7 > 0$ and $j_5 = j_5(\Omega) \in \mathbb{N}$ such that for all $\theta \in (0, \pi/2)$ corresponding to a $j \geq j_5$ reflection orbit emanating from $(x_0, \xi(\theta))$ which makes approximately one rotation, we have $\d x_j / \d \theta \geq C_7 j$. Similarly, for $\theta \in (-\pi/2, 0)$ corresponding to an orbit of $j$ reflections and approxmately one rotation, we have $\d x_j / \d \theta \leq - C_7 j$.
\end{lemm}
\begin{proof}
Denote $x_j = \pi_1 \beta^j(x_1,\alpha_1)$, where $\pi_1$ is the projection onto the first component. Then $x_j$ depends only on $(x_0, \xi)$, or equivalently $(x_0, \theta)$ and can be written as the composition
\begin{align*}
x_j = \pi \beta^j(x_1(x_0,\theta), \alpha_1(x_0,\theta)).
\end{align*}
As $\theta$ increases or decreases from $0$, we have
\begin{align}\label{chain rule}
\frac{\d x_j}{\d \theta} = \frac{\d x_j}{\d x_1} \frac{\d x_1}{\d \theta} + \frac{\d x_j}{\d \alpha} \frac{\d \alpha}{\d \theta}.
\end{align}
Recall from Lemma \ref{boundary monotonicity} that
\begin{align}\label{large derivative}
C_3 j \leq \frac{\d x_j}{\d \alpha} \leq C_4 j,
\end{align}
for $j$ large and positive constants $C_3$ and $C_4$. This derivative can be made arbitrarily large by choosing $j$ accordingly. Also, Lemma \ref{L3} showed that all orbits making approximately one rotation with $j$ reflections at the boundary have angles of reflection in the range $C_1/j \leq \alpha \leq C_2/j$ for positive constants $C_1$ and $C_2$. We can now use Lemma \ref{L1}, which showed that $|\d \alpha/\d \theta| \geq c > 0$ independently of $j$ for all $\theta$ producing an orbit of approximately one rotation in $j$ reflections. Using Lazutkin coordinates, we also have that
\begin{align}
\frac{\d x_j}{\d x_1} \frac{\d x_1}{\d \theta} = (1 + O(1/j)) \frac{\d x_1}{\d \theta}.
\end{align}
In graph coordinates (see proof of Lemma \ref{L1}), one can calculate that
\begin{align}
\begin{split}
\left|\frac{\d x_1}{\d \theta}\right| &= \left| \frac{\d x_1}{\d x} \frac{\d x}{\d \theta}\right| = \sqrt{1 + (f'(x))^2} \frac{x f'(x) + b - f}{x^2 + (b - f(x))^2}\\
& \leq 1 + O(b) \leq 2
\end{split}
\end{align}
which is bounded independently of $j$. Hence, the term \eqref{large derivative} dominates in the expression \eqref{chain rule} and $x_j$ winds monotonically around $\d \Omega$.
\end{proof}

\noindent At the endpoints $\theta = \pm \pi/2$, it is clear that the angle of reflection at $x_{\pm 1}$ is $\pi/2$ since the distance curves are parallel. Angles outside the range $\theta \in (-\pi/2, \pi/2)$ correspond to clockwise orbits and Lemma \ref{L3} showed that for an orbit making approximately one rotation and $j$ reflections at the boundary, all angles of reflection satisfy $C_1/j \leq \alpha_k \leq C_2/j$, $(1 \leq k \leq j)$. As $j$ reflection orbits emanating from $\xi_0^+$ (corresponding to $\theta = 0$) make only a quarter rotation, it is clear that the collection of covectors at $x_0$ whose trajectories make approximately one rotation in $j$ reflections consists of two connected components in $S_{x_0}^*\Omega$.
\begin{def1}
The positive admissible cone $C_{x_0, +}^*(\Omega; j)$ at $x_0$ is defined to be the set of homogeneous extensions to $T_{x_0}^* \Omega$ of the two components in $S_{x_0}^*(\Omega)$ described above. The negative admissible cone $C_{x_0, -}^*(\Omega; j)$ at $x_0$ is defined by the same property for orbits making $j$ reflections in the clockwise direction. Their union is denoted by $C_{x_0}^*(\Omega;j) = C_{x_0, +}^*(\Omega;j) \cup C_{x_0, -}^*(\Omega;j)$. See Figure \ref{picture of counterclockwise perturbation} for an illustration of $C_{x_0}^*(\Omega;j)$.
\end{def1}

\subsection{Intersection points}
We are finally ready to show that the intersection points of the last link with the distance curve of $y$ wind around monotonically as we twist $\xi$ in either direction. We first explain why the last link necessarily intersects $d^{-1}(\text{dist}(y, \d \Omega) )$ twice.
\begin{lemm}\label{2int}
There exists $j_6 = j_6(\Omega) \in \mathbb{N}$ such that for $j \geq j_6$ and any $N \geq 4$, the distance curve $\{x \in \Omega: \text{dist}(x,\d \Omega) = j^{-N}\}$ is intersected exactly twice by any link emanating from the boundary at an angle greater than or equal to $C_1/j$, with $C_1$ the constant appearing in Lemma \ref{L3}.
\end{lemm}
\begin{proof}
For each point $p$ in the distance curve $\{z \in \Omega : d(z, \d \Omega) = j^{-N} \}$, the tangent line $T_p \{z \in \Omega : d(z, \d \Omega) = j^{-N} \}$ intersects $\d \Omega$ exactly twice by convexity. Recall that Lemma \ref{L3} gave $C_1/j \leq \alpha_k \leq C_2/ j$ for each $1 \leq k \leq j$. We now show that the angles of reflection for links which are tangent to the distance curve $d^{-1}(j^{-N})$ are of order $O(j^{-N/2})$. In graph coordinates, $\d \Omega$ is locally given by $(x,f(x))$ where $f(x) = \kappa x^2 + R(x)$, $\kappa$ is the curvature of $\d \Omega$ at the point corresponding to $(0,0)$ and $R(x) = O(x^3)$. By rescaling, we may assume that $\kappa = 1$. The distance curve $d^{-1}(j^{-N})$ can be parametrized in graph coordinates by
\begin{align*}
d^{-1}(j^{-N}) = \left\{ (t, f(t)) + r N(t): t \in \mathbb{R}  \right\},
\end{align*}
where $r = j^{-N}$ and
\begin{align*}
N(t) = \frac{(-f'(t), 1)}{\sqrt{1+ (f'(t))^2}}
\end{align*}
is the unit normal to the graph at $(t,f(t))$. There exist precisely two lines through the origin which are tangent to $d^{-1}(j^{-N})$ in graph coordinates. The positive parameter $t$ corresponding to a point of tangency satisfies
\begin{align}\label{slopes}
\frac{f(t) + r/(1 + (f'(t) )^2)^{1/2}}{t - rf'(t)/(1 + (f'(t) )^2 )^{1/2} } = \frac{\d_t (f(t) + r/(1 + (f'(t) )^2 )^{1/2} ) }{\d_t (t - rf'(t) /(1 + (f'(t) )^2)^{1/2})},
\end{align}
where the lefthand side of \eqref{slopes} is the slope of a line connecting the origin to a point on $d^{-1}(r)$ and the righthand side is the slope of the tangent to $d^{-1}(r)$. Solving for $t$ in terms of $r$ in equation \eqref{slopes} yields $f'(t) t - f(t) = t^2 + O(t^3) = O(r)$. Plugging this in to the righthand side of \eqref{slopes}, we see that the angle of tangency satisfies $\tan (\alpha_{\text{tangency}}) = O(r^{1/2})$ and hence $\alpha_{\text{tangency}} = O(j^{-2})$ if $N \geq 4$. The proof is complete by noting that the angles $\alpha_k$ coming from orbits making approximately one rotation and $j$ reflections are bounded below by $C_1/j > O(j^{-4/2})$ for $j$ sufficiently large.

\end{proof}

\begin{lemm}\label{clockwise int points}
Under the hypotheses of Lemma \ref{2int}, denote by $w_1$ and $w_2$ the arclength coordinates on $d^{-1}(d(y))$ corresponding to the intersection points of $d^{-1}(d(y))$ and the last link of the $j$ reflection orbit emanating from $x$. There exist $C_8 > 0$ and $j_7 = j_7(\Omega) \in \mathbb{N}$ such that if $j \geq j_7$ and $x,y$ are $O(1/j^4)$ close the diagonal of the boundary, then $|\d w_i/ \d\theta| \geq C_8 j$ ($i = 1,2$) for all $\theta$ corresponding to $\xi(\theta)$ in the cone of admissible covectors at $x_0$.
\end{lemm}
\begin{proof}
In graph coordinates, $\d \Omega$ is again locally parametrized by $(t,f(t))$, where $f(t) = \kappa t^2 + R(t)$, $\kappa$ is the curvature of $\d \Omega$ at the point corresponding to $(0,0)$ and $R(t) = O(t^3)$. By rescaling, we may assume that $\kappa = 1$. The distance curve $d^{-1}(d(y))$ on which $y$ lies can be locally parametrized in graph coordinates by
\begin{align}\label{d parametrization}
d^{-1}(d(y)) = \left\{ (t, f(t)) + r N(t)  \right\},
\end{align}
where $r = \text{dist}(y,\d \Omega) = O(j^{-4})$ and
\begin{align*}
N(t) = \frac{(-f'(t), 1)}{\sqrt{1+ (f'(t))^2}}
\end{align*}
is the unit normal to the graph at $(t,f(t))$. The final link of the orbit connecting $x_j(\theta) \in \d \Omega$ to $x_{j+1}(\theta) \in \d \Omega$ can be parametrized by the line
\begin{align}\label{last link line}
Y = \tan (\alpha_j(\theta) + \beta_j(\theta)) X + f(t_j(\theta)) - \tan(\alpha_j(\theta) + \beta_j(\theta)) t_j(\theta),
\end{align}
where
\begin{align*}
\beta_j = \arccos \left( \frac{1}{\sqrt{1 + (f'(x_j))^2}}  \right)
\end{align*}
is the angle of the tangent to the graph with the horizontal axis and $t_j = t_j(\theta) \in \mathbb{R}$ is defined implicitly by the equation $(t_j,f(t_j)) = x_j \in \d \Omega$. Noting that the Cartesian coordinates of the parametrization of $d^{-1}(d(y))$ satisfy
\begin{align}
\begin{split}
X &= t - \frac{rf'(t)}{\sqrt{1+ (f'(t))^2}},\\
Y &= f(t) + \frac{r}{\sqrt{1+ (f'(t))^2}},
\end{split}
\end{align}
we plug these into equation \eqref{last link line} to obtain
\begin{align}\label{implicit eqn}
\begin{split}
\left(f(t) + \frac{r}{\sqrt{1+ (f'(t))^2}}\right) &= \tan (\alpha_j(\theta) + \beta_j(\theta)) \left(t - \frac{rf'(t)}{\sqrt{1+ (f'(t))^2}} \right)\\
&+ f(t_j(\theta)) - \tan(\alpha_j(\theta) + \beta_j(\theta)) t_j(\theta).
\end{split}
\end{align}
\noindent By Lemma \ref{2int}, there exist precisely two solutions $t = z_1(\theta), z_2(\theta)$ of equation \eqref{implicit eqn} which correspond to the two intersection points of the last link with the distance curve $d^{-1}(d(y))$. Plugging $z_{i}(\theta)$ into equation \eqref{implicit eqn}, differentiating in $\theta$ and evaluating at $\theta_0$ corresponding to $x_j(\xi(\theta_0))$ at the origin, we obtain
\begin{align}\label{zprime}
\begin{split}
f'(z_i)& z_i' - \frac{r f'(z_i) f''(z_i) z_i'}{(1 + f(z_i)^2)^{3/2}} =   (\alpha_j' + \beta_j') (1 + \tan^2 \alpha_j) \left(z_i - \frac{rf'(z_i)}{\sqrt{1+ (f'(z_i))^2}} \right)\\
&+  \tan \alpha_j \left(z_i' - \frac{rf''(z_i) z_i'}{\sqrt{1+ (f'(z_i))^2}}  + \frac{r(f'(z_i))^2 f''(z_i) z_i'}{(1 + (f'(z_i))^2 )^{3/2}} \right) + 0 - \tan \alpha_j t_j'.
\end{split}
\end{align}
We now collect all terms with $z_i'$ in \eqref{zprime}:
\begin{align}\label{divide through}
\begin{split}
A_i z_i' &=   (\alpha_j' + \beta_j') (1 + \tan^2\alpha_j ) \left(z_i - \frac{rf'(z_i)}{\sqrt{1+ (f'(z_i))^2}} \right) - \tan \alpha_j t_j',
\end{split}
\end{align}
where
\begin{align*}
A_i = f'(z_i)  - \frac{r f'(z_i) f''(z_i)}{(1 + f(z_i)^2)^{3/2}} -  \tan \alpha_j \left(1 - \frac{rf''(z_i)}{\sqrt{1+ (f'(z_i))^2}}  + \frac{r(f'(z_i))^2 f''(z_i) }{(1 + (f'(z_i))^2 )^{3/2}}\right).
\end{align*}
Before dividing through equation \eqref{divide through} by $A_i$, we Taylor expand $A_i$ in powers of $1/j$ to show that it is nonvanishing:
\begin{align}
A_i = 2 z_i - \tan \alpha_j + O(j^{-2}).
\end{align}
At this point, it is important to distinguish between $z_1$ and $z_2$. If $z_1 < z_2$, then insisting that $y$ be $O(j^{-4})$ close to the boundary amounts to setting $z_1 = O(j^{-4})$ and $z_2 = t_{j+1} + O(j^{-4})$. Here, $t_{j+1} = t_{j+1}(\theta)$ is defined implicitly by the equation $(t_{j+1}, f(t_{j+1})) = x_{j+1} \in \d \Omega$. When $i = 1$, $A_1 = - \alpha_j + O(j^{-4}) \leq -C_1 / (2j)$ is nonvanishing for $j$ sufficiently large. For $i = 2$, note that $\tan \alpha_j = f(z_2)/z_2 + O(j^{-4})$, which implies that $A_2 = z_2 + O(j^{-2})$ is again nonvanishing. Hence, we may divide through equation \eqref{divide through} by $A_i$. To estimate the righthand side of \eqref{divide through}, we calculate that
\begin{align}
\beta_j' = \frac{f'(t_j)}{\sqrt{1 - \frac{1}{1 + (f'(t_j))^2 }}} \frac{ f''(t_j) t_j '}{(1 + (f(t_j)^2))^{3/2}}.
\end{align}
The first factor above is well defined by continuity at $t_j = 0$ and equals $1$. Hence, $\beta_j' = 2 t_j' + O(1/j)$. If $i = 1$, equation \eqref{divide through} becomes
\begin{align}
(- \tan \alpha_j  +O(j^{-4})) z_1' = - \tan \alpha_j t_j' + O(j^{-3}),
\end{align}
which implies that $z_1' = t_j' + O(j^{-2})$. When $i = 2$, equation \eqref{divide through} becomes
\begin{align}
(z_2 + O(j^{-2}))z_2' = t_j' z_2 + O (1/j),
\end{align}
which again implies that $z_2' = t_j' + O(1/j)$. Recall that Lemma \ref{winding in xi} showed $x_j'$ is of order $j$ in Lazutkin coordinates when $\xi(\theta) \in C_{x_0}^*(\Omega;j)$. To compare Lazutkin coordinates with $t_j$, note that in graph coordinates, the arclength parameter $ds$ on $\d \Omega$ is given by
$$
ds = (1 + (f'(t))^2)^{1/2} dt.
$$
As arclength coordinates and Lazutkin coordinates are comparable independently of $j$, we conclude that $t_j'$ and hence $z_i'$ are also of order $j$. We conclude the proof by noting that the arclength parameter $ds'$ on $d^{-1}(d(y))$ in the parametrization \eqref{d parametrization} is given by
$$
ds' = \left((1 + (f'(t))^2)^{1/2} + O(r)\right) dt,
$$
which is also comparable to $dt$ independently of $j$.
\end{proof}

\subsection{Clockwise orbits}
In Lemma \ref{clockwise int points}, it was shown that there are precisely $4$ counterclockwise orbits connecting $x = x_0$ to $y$ in $j$ reflections and approximately one rotation. The only constraint on $x$ and $y$ was that they were confined to a $O(j^{-4})$ neighborhood of the diagonal of the boundary. By reflecting $\Omega$ through the verticle axis, one obtains another smooth strictly convex domain and the reflections of $x$ and $y$ remain $O(j^{-4})$ close to the diagonal of the boundary. Hence, the same procedure produces exactly $4$ counterclockwise orbits of $j$ reflections from $x$ to $y$ making approximately one rotation in the reflected domain. Reflecting back through the vertical axis carries these $4$ orbits to clockwise orbits in the original domain. This completes the proof of Theorem \ref{8 orbit lemma}.

\begin{rema}\label{Uj}
We can take $j_0 = j_0(\Omega) = \max\{j_i: 2 \leq i \leq 7\}$, with $j_i$ as they appear in Lemmas \ref{L3}, \ref{L1}, \ref{boundary monotonicity}, \ref{winding in xi}, \ref{2int}, and \ref{clockwise int points}. Similarly, we can choose a uniform constant $C_0$ as in the statement of Lemma \ref{8 orbit lemma}. The tubular neighborhood $U_j$ referenced in Theorem \ref{HRP} can be taken to be $\{(x,y) \in \Omega \times \Omega: \text{dist}((x,y), \Delta \d \Omega) \leq j^{-4}\}$ for $j \geq j_0$, where $\Delta: \d \Omega \times \d \Omega$ is the diagonal embedding $x \mapsto (x,x)$. For $(x,y) \in U_j$, both $x$ and $y$ satisfy the conditions of Lemma \ref{clockwise int points}. The cone bundle $C_{x}^*(\Omega;j)$ is well efined whenever $\text{dist}(x,\d \Omega) = O(j^{-4})$.
\end{rema}

\begin{rema}
The proof of Theorem \ref{8 orbit lemma} in this section could actually be extended to a larger region of validity. In particular, the same methods allow us to prove the existence of $8$ orbits of rotation number $k/j$ for $j$ sufficiently large in terms of $k$, connecting points in a comparable open neighborhood of the diagonal of the boundary. Additionally, we could also allow $x$ and $y$ to be further away from the diagonal as long as they are both sufficiently close to the boundary. However, we are only concerned with near diagonal terms in this paper for the purposes of deriving trace formulas.
\end{rema}

\begin{rema}\label{boundary points}
If $x \in \d \Omega$, then $\xi_0^\pm$ are tangent to the boundary and perturbing $\xi$ away from $\xi_0^+$ in the clockwise direction is no longer well defined. Rather, it is equivalent to reflecting and then rotating $\xi$ in the counterclockwise direction. Similarly, $\xi_0^-$ cannot be rotated in the counter clockwise direction. Each of these restrictions reduces the number of $j$ reflection orbits to $y$ in appoximately one rotation by $2$. If additionally $y \in \d \Omega$, the final link only makes one intersection with $d^{-1}(0) = \d \Omega$ so there are only two orbits of $j$ reflections connecting $x$ to $y$ in approximately one rotation. One is in the counterclockwise direction while the other is in the clockwise direction. In particular, when $x = y$, there is a \textit{unique} geodesic loop (up to parametrization) of $j$ reflections and exactly one rotation (i.e. rotation number $1/j$). The proof of these claims is much simpler, as it only requires Lemma \ref{boundary monotonicity}.
\end{rema}

\section{A parametrix for the wave propagator}\label{A parametrix for the wave propagator}
\noindent In this section, we use microlocal analysis to obtain an oscillatory integral parametrix for the wave propagator in an open subset of $\mathbb{R} \times \Omega \times \Omega$ which contains the diagonal $\mathbb{R} \times \Delta(\d \Omega)$. The wave kernel $e^{it\sqrt{- \Delta}}$ is actually \textit{not} a Fourier integral operator (FIO) near the tangential rays (see \cite{AndersonMelrose}, \cite{MelroseTaylor}), so we microlocalize the wave kernels near periodic transversal reflecting rays. We begin by reviewing FIOs and Chazarain's parametrix for the wave propagator in Sections \ref{Fourier Integral Operators} and \ref{Chazarain's Parametrix}. In Section \ref{Oscillatory integral representation}, we then cook up oscillatory integrals for each term in Chazarain's parametrix, which microlocally approximate the wave propagator near the orbits described in Theorem \ref{8 orbit lemma}. This approach is inspired by the formulas in \cite{MaMe82} and will rely on the symbol calculus in Section \ref{Chazarain's Parametrix}.

\subsection{Fourier Integral Operators}\label{Fourier Integral Operators}
Let $X$ and $Y$ be open sets in $\mathbb{R}^{n_X}$ and $\mathbb{R}^{n_Y}$ respectively. If $a \in S^\mu_{1, 0}(X \times \mathbb{R}^N)$ is a classical symbol of order $\mu$ and $\Theta \in C^\infty(X  \times \mathbb{R}^N)$ is a nondegenerate phase function, then the linear operator
$$
A(u) = \int_X \int_{\mathbb{R}^N} e^{i \Theta(x,\theta)} a(x,\theta) u(x) \,d\theta dx
$$
is called a Lagrangian or Fourier integral distribution on $X$. Recall that a continuous linear operator $A : C_0^\infty(Y) \to \mathcal{D}'(X)$ has an associated Schwartz kernel $K_A \in \mathcal{D}'(X \times Y)$. If $K_A$ is given by a locally finite sum of Lagrangian distributions on $X \times Y$, then we say $A$ is a Fourier integral operator (FIO). One can then show that the wavefront set of the kernel is contained in the image of the map $\iota_{\Theta}: (x, y) \mapsto (x, y, d_x \Theta, d_y\Theta)$ when restricted to the critical set $C_\Theta := \{d_\theta \Theta = 0\}$. The image of $\iota_\Theta$ is in fact a conic Lagrangian submanifold $\Lambda_\Theta \subset T^* (X \times Y)$ and the map $\iota_{\Theta}$ is a local diffeomorphism from $C_\Theta$ onto $\Lambda_\Theta$. In this case, we say ``$\Theta$ parametrizes $\Lambda_\Theta$.'' The canonical relation or wavefront relation of $A$ is defined by
$$
WF'(A) = \{(x,\xi), (y, \eta) : (x,y,\xi, -\eta) \in WF(K_A)\} \subset T^*X \times T^*Y,
$$
and describes how the operator $A$ propagates singularities of distributions on which it acts. If $\omega_X$ and $\omega_Y$ denote the natural symplectic forms on $T^*X$ and $T^*Y$ respectively, one can more invariantly consider FIOs associated to general conic Lagranigan submanifolds $\Lambda \subset T^*X \times T^*Y$ (canonical relations), with respect to the symplectic form $\omega_X - \omega_Y$. The notion of a principal symbol for Fourier integral operators is more subtle than that for pseudodifferential operators. The principal symbol of $A$ is a half density on $\Lambda_\Theta$ given in terms of the parametrization $\iota_\Theta$:
\begin{align}\label{principal symbol of an FIO}
e = {\iota_\Theta}_* (a_0 |dC_\Theta|^{1/2}),
\end{align}
where $a_0$ is the leading order term in the asymptotic expansion for $a$ and $|dC_\Theta|^{1/2}$ is the half density associated with the Gelfand-Leray form on the level set $\{d_\theta \Theta = 0\}$. Here, we have ignored Maslov factors coming from the Keller-Maslov line bundle over $\Lambda_\Theta$. These are nonzero factors $e^{i\sigma \pi/4}$ ($\sigma$ is known as the Maslov index) which appear in front of the principal symbol as a result of the multiplicity of phase functions parametrizing the canonical relation $\Lambda_\Theta$, possibly in different coordinate systems. While these factors allow the principal symbol to be defined in a more geometrically invariant way, we defer computation of the Maslov indices until Section \ref{Maslov factors}. For a more thorough reference on the global theory of Lagrangian distributions, see \cite{Du96}, \cite{Ho71}, \cite{DuHo72}, and \cite{Ho485}. The order of a Fourier integral operator is defined in such a way that when two Fourier integral operators' canonical relations intersect transversally, then the composition is again a Fourier integral operator and order of the composition is the sum of the orders:
\begin{align}\label{order of a lagrangian distribution}
\text{order}(A) = m = \mu + \frac{1}{2}N - \frac{1}{4}(n_X + n_Y).
\end{align}
Recall that here, $n_X$ and $n_Y$ are the dimensions of $X$ and $Y$ respectively. In this case, we write $A \in I^m(X \times Y, \Lambda)$. This convention on orders also generalizes that of pseudodifferential operators, where $X = Y$ and $m = \mu$ coincides with the order of the corresponding symbol class. Sufficient conditions which guarantee that the composition exists are clean or transversal intersection of the two operators' canonical relations. In general, composition of Fourier integral operators and the associated symbol calculus is somewhat complicated, but we will not directly use the composition formula in what follows.

\subsection{Chazarain's parametrix}\label{Chazarain's Parametrix}
The parametrix developed by Chazarain in \cite{Ch76} provides a microlocal description of the wave kernels near periodic transversal reflecting rays. The parametrices for
$$
E(t) = \cos t \sqrt{- \Delta}, \qquad \qquad S(t) = \frac{\sin t \sqrt{- \Delta}}{\sqrt{-\Delta}}
$$
are constructed in the ambient Euclidean space $\mathbb{R} \times \mathbb{R}^n \times \mathbb{R}^n$. We only consider $S(t)$, as the formula for $E(t)$ is easily obtained from that of $S(t)$ by differentiating in $t$. We write $E(t,x,y)$ and $S(t,x,y)$ for the Schwartz kernels of $E(t)$ and $S(t)$ respectively. Following the work in \cite{Ch76}, we can find a Lagrangian distribution
\begin{align}\label{Chazarain sum}
\widetilde{S}(t,x,y) = \sum_{j = - \infty}^{\infty} S_j(t,x,y), \qquad S_j \in I^{-5/4}(\mathbb{R} \times \mathbb{R}^n \times \mathbb{R}^n, \Gamma_\pm^j),
\end{align}
which approximates $S(t)$ microlocally away from the tangential rays modulo a smooth kernel. We will describe the canonical relations $\Gamma_\pm^j$ momentarily and in particular, show that the sum in \eqref{Chazarain sum} is locally finite. We first explain what is meant by approximating $S(t)$ ``microlocally away from the tangential rays.'' In general, two distributions $f,g \in \mathcal{D}'(\mathbb{R}^n)$ are said to agree microlocally near a closed cone $\Lambda_1 \subset T^*\mathbb{R}^n$ if $WF(u - v) \cap \Lambda_1 = \emptyset$. Similarly, using the language from Section \ref{Fourier Integral Operators}, two operators $A, B : C^\infty(Y) \to C^\infty(X)$ are said to agree microlocally near a closed cone $\Lambda_2 \subset T^*X \times T^*Y$ if $WF'(A-B) \cap \Lambda_2 = \emptyset$. This second notion is what we will use to say that our parametrix approximates $S(t)$ microlocally near the canonical relations $\Gamma_{\pm}^j$.
\\
\\
To describe the canonical relations precisely, we first introduce some notation following the presentations in \cite{Ch76} and \cite{HeZe12}. As the Euclidean wave operator $\Box_{\mathbb{R}^2}$ factors into $(\d_t - \sqrt{\Delta})(\d_t + \sqrt{-\Delta})$, there are two Hamiltonians corresponding to the symbol $\pm |\eta|$ of $\pm \sqrt{-\Delta}$. Let $H_\pm(y,\eta) = \pm |\eta|$ and $g^t$ be the Hamiltonian flow, i.e. the flow map associated to the system
\begin{align*}
\begin{cases}
\d_t y = \frac{\d H_\pm}{\d \eta},\\
\d_t \eta = -\frac{\d H_\pm}{\d y},
\end{cases}
\end{align*}
which is in fact just the reparametrized geodesic flow on $\mathbb{R}^2$. For $(y,\eta) \in T^*\Omega$ (or $T_{\d \Omega}^*\mathbb{R}^2$ such that $\eta$ is transversal to the boundary and inward pointing), recall that in Section \ref{Billiards} we defined
\begin{align*}
\begin{split}
t_\pm^1 &= \inf_{t > 0} \{t: \pi_1 g^{t}(y, \eta) \in \d \Omega \},\\
t_\pm^{-1} &= \sup_{t < 0} \{t: \pi_1 g^{t}(y,\eta) \in \d \Omega \},
\end{split}
\end{align*}
where $\pi_1$ is projection onto the spatial variable. Note that $t_\pm^{\pm 1}$ depend on $(y,\eta)$. We also have $t_\mp^{1} = -t_\pm^{-1}$. We then set
\begin{align*}
\begin{split}
\lambda_\pm^{1}(y, \eta) &= g^{\pm t_\pm^1}(y,\eta),\\
\lambda_\pm^{-1}(y, \eta) &= g^{\pm t_\pm^{-1}}(y,\eta).
\end{split}
\end{align*}
Also define $\widehat{\lambda_\pm^{1}(y,\eta)}$ to be the reflection of $\lambda_\pm^{1}(y,\eta)$ through the cotangent line at the boundary. In other words, $\widehat{\lambda_\pm^{1}(y,\eta)}$ and $\lambda_\pm^{1}(y,\eta)$ have the same cotangential components but opposite conormal components so that $\widehat{\lambda_\pm^{1}(y,\eta)}$ is inward pointing. The point $\widehat{\lambda_\pm^{-1}(y,\eta)}$ is defined analogously. We can then inductively define $t_\pm^j$ and $t_\pm^{-j}$ by the formulas
\begin{align*}
\begin{split}
t_\pm^j &= \inf_{t > 0} \{t : \pi_1 g^{\pm t}(\widehat{\lambda_\pm^{j-1}(y,\eta)}) \in \d \Omega \},\\
t_\pm^{-j} &= \sup_{t < 0} \{t : \pi_1 g^{\pm t}(\widehat{\lambda_\pm^{-(j-1)}(y,\eta)}) \in \d \Omega \}.
\end{split}
\end{align*}
The total travel time after $j$ reflections is defined by
\begin{align*}
T_\pm^j = \begin{cases}
\sum_{k = 1}^j t_\pm^k & j > 0\\
0 & j = 0\\
\sum_{k = j}^{-1} t_\pm^k & j < 0.
\end{cases}
\end{align*}
To study how the fundamental solution of $\Box_{\Omega}$ behaves at $\d \Omega$ when we impose boundary condtions, we propagate the intial data by the free wave propagator on $\mathbb{R}^2$, restrict it to the boundary, reflect, and then propagate again.  If such a construction is continued for $j \in \mathbb{Z}$ reflections at the boundary, it is shown in \cite{Ch76} that the FIOs $S_j$ must have canonical relations
$$
\Gamma_\pm^j = \begin{cases}
(t, \tau, g^{\pm t}(y,\eta), y, \eta): \tau = \pm |\eta| & j =0,\\
(t, \tau, g^{\pm (t - T_\pm^j(y,\eta))}(\widehat{\lambda_\pm^j(y,\eta)}), y, \eta)): \tau = \pm |\eta| & j \in \mathbb{Z} \backslash \{0\}.
\end{cases}
$$
\noindent Since $\widetilde{S}(t)$ is a microlocal parametrix, the canonical relation of the true solution operator $S(t)$ is also
$$
\Gamma = \bigcup_{j \in \mathbb{Z}, \pm} \Gamma_{\pm}^j.
$$
\noindent Here, $j > 0$ and $j < 0$  correspond to reflections on the inside and outside of the boundary. In fact, there are four modes of propagation associated to the canonical relations $\Gamma_\pm^j$, corresponding to $\pm \tau \geq 0$ (forwards and backwards time) and $\pm j \geq 0$ (inside and outside bounces). See \cite{HeZe12} for an explicit local model in the case of one reflection, where $\Omega$ is replaced by a half plane. We see that a point $(t, \tau, x, \xi, y, \eta)$ belongs to the canonical relation $\Gamma_{\pm}^j \subset T^*(\mathbb{R} \times \mathbb{R}^2 \times \mathbb{R}^2)$ if the broken geodesic of $j$ reflections emanating from $(y,\eta)$ passes through $(x,\xi)$. For a more thorough discussion of Chazarain's parametrix, see \cite{Ch76} and \cite{GuMe79b}.
\\
\\
\noindent To be precise, Chazarain actually showed that there exists FIOs $S_j$ such that the sum in \eqref{Chazarain sum} is a parametrix for the wave propagator $S(t)$ with canonical relation
$$
\Gamma = \bigcup_{j \in \mathbb{Z}, \pm} \Gamma_{\pm}^j.
$$
However, the canonical relations were never parametrized by explicit phase functions and the principal symbols of the operators $S_j$ were not computed in \cite{Ch76}. For the remainder of this section, we concern ourselves with the task of explicitly computing them in terms of geometric and dynamical data associted to the billiard map.

\subsection{Oscillatory integral representation}\label{Oscillatory integral representation}
In this section, we cook up an oscillatory integral such that microlocally near the canonical relations $\Gamma_\pm^j$,
$$
S_j(t,x,y) = \int_{-\infty}^\infty e^{i\Theta_j(t,\tau, x,y)} {a_j}(\tau, x,y)\,d\tau,
$$
where $S_j$ is the $j^{th}$ term in Chazarain's parametrix corresponding to a wave with $j$ reflections and $a_j \in S_{\text{cl}}^\mu$ is a classical symbol of order $\mu$. We will only compute the principal symbol and use $L.O.T$ to denote lower order terms in the sense of Lagrangian distributions in what follows. Due to the presence of different Maslov factors on each branch of $\Gamma_{\pm}^j$ corresponding to $\pm \tau > 0$ (cf. Sections \ref{Fourier Integral Operators} and \ref{Maslov factors}), it is actually more convenient to find operators
\begin{align}\label{osc int}
S_{j, \pm}(t,x,y) = \int_{0}^\infty e^{i\Theta_{j,\pm}(t,\tau, x,y)} { a_{j,\pm}}(\tau, x,y)\,d\tau,
\end{align}
so that $S_j = S_{j,+} + S_{j,-}$ and the phase functions associated to $S_{j,\pm}$ paramaterize $\Gamma_{+}^j$ and $\Gamma_-^j$ individually.
\\
\\
We first make precise the notion of microlocalized FIOs. We would like to microlocalize $S(t)$ near periodic orbits of rotation number $1/j$. Oftentimes, it is required that such orbits be nondegenerate, in the sense that $1$ is not an eigenvalue of the linearized Poincar\'e map. However, this assumption is not needed for our trace formulas, which work both for simple nondegenerate orbits as well as degenerate orbits coming in one parameter families as in the case of an ellipse (cf. Theorem \ref{Rational caustic theorem} and Corollary \ref{Ellipse corollary}). We recall that $\Omega$ is said to satisfy the concoincidence condition \eqref{NCC} if there exist no periodic orbits of rotation number $m/n$, $m \geq 2$ in a sufficiently small neighborhood of $|\d \Omega|$. In this case, for $j$ sufficiently large, periodic orbits of rotation number $1/j$ come in isolated families. This follows from the results in \cite{MaMe82}, which show that for $j$ sufficiently large, no two orbits of distinct rotation numbers $1/j$ and $1/k$ can have the same length. It can be shown that periodic billiard orbits arise as critical points of the loop function $\Psi_j(q,q)$ (see \cite{Vig18}, \cite{HeZe19}). As in the statement of Theorem \ref{Main theorem}, denote $t_j = \inf_{q \in \d \Omega}$ and $T_j = \sup_{q \in \d \Omega}$. If $\Omega$ satisfies the noncoincidence condition \ref{NCC}, we can find a smooth cutoff function $\chi_1(t)$ which is identically equal to $1$ on an open nieghborhood of $[t_j,T_j]$ and vanishes in a neighborhood of all other $L \in \text{Lsp}(\Omega)$. As noted in Section \ref{Chazarain's Parametrix}, each propagator $S_j$ has canonical relations $\Gamma_{\pm}^j$. Denote by $\chi_2$ a smooth cutoff function which is identically equal to $1$ on $\Gamma^j = \cup_\pm \Gamma_\pm^j$, vanishing near the gliding rays $T^* \d \Omega$ and conic in the fiber variables $\tau, \xi$ and $\eta$. Quantizing $\chi_2$ gives a pseudodifferential operator which microlocalizes near the support of $\chi_2$. For a reference, see Chapter 18 of \cite{Ho385}. We call such an operator a microlocal cutoff on $\Gamma_\pm^j$. The composition $\chi_1(t)\chi_2(t,x,y, D_t, D_x, D_y) S(t)$ is then smoothing away from the geodesic loops of rotation number $1/j$ and the trace of this composition is equal to the wave trace modulo $C^\infty$ in an open neighborhood of $[t_j,T_j]$.
\\
\\
We now use Theorem \ref{8 orbit lemma} to find suitable phase functions $\Theta_{j,\pm}$ which parametrize $\Gamma_\pm^j$. Define phase functions $\Theta_j^k$ by the formula
$$
\Theta_{j,\pm}^k(t,\tau, x,y) = \pm \tau(t - \Psi_j^k(x,y)),
$$
where $\Psi_j^k$ are given in Definition \ref{psi jk}. We then have,

\begin{lemm}\label{parameterization of cannonical relations}
The phase functions $\Theta_{j,\pm}^k(t,\tau,x,y)$ are smooth in an open neighborhood of the diagonal of the boundary and locally parametrize the canonical relations $\Gamma_\pm^j$. In particular, the fibers of both $\Gamma_+^j$ and $\Gamma_-^j$ lying over this neighborhood are unions of $8$ connected components, which we denote by $\Gamma_\pm^{j,k}$ corresponding to $1 \leq k \leq 8$ as in Definition \eqref{psi jk}.
\end{lemm}
\noindent A proof of Lemma \ref{parameterization of cannonical relations} can be found in \cite{Vig18} and is based on the fact that the length functional is a generating function for the billiard map. For a fixed $j$, the multiple components $\Gamma_\pm^{j,k}$ of the canonical relation indexed by $k$ can be seen to intersect over the boundary. As our parametrix for the propagators described in \cite{Ch76} are in fact modified by microlocal cutoffs supported away from the tangential rays in $S^*\d \Omega$, we may assume $\Gamma_{\pm}^{j,k}$ are smooth nonintersecting Lagrangian submanifolds. It would be interesting to study mapping properties of the operators $S_j^k$ as the canonical relations degenerate near the glancing set in future work. We now want to derive an explicit oscillatory integral representation \eqref{osc int} for $S(t)$ in a specific coordinate system adapted to $\Omega$. We first need to better understand the forwards and backwards symbols on $\Gamma$.
\begin{prop}\label{symbol prop}
Let $e_\pm$ denote the principal symbol of ${S}(t)$ on $\Gamma = \bigcup_{j \in \mathbb{Z}, \pm} \Gamma_\pm^j$. Then, we have
\begin{align*}
e_\pm &= \frac{(-1)^j}{2 \tau i} |dt \wedge dy \wedge d\eta|^{1/2} \qquad \text{(Dirichlet)}\\
e_\pm &= \frac{1}{2 \tau i} |dt \wedge dy \wedge d\eta|^{1/2} \qquad \text{(Neumann/Robin)}
\end{align*}
where $|dt \wedge dy \wedge d\eta|^{1/2}$ is the canonical half density.
\end{prop}
\begin{proof}
For a proof in the Dirichlet/Neumann case, see \cite{HeZe12}. For the Robin boundary conditions, see \cite{Vig18}.
\end{proof}

\noindent As Proposition \ref{symbol prop} gives $e_\pm = \frac{1}{2 \tau i} |dt \wedge dy \wedge d\eta|^{1/2}$ and we now know that the phase functions $\Theta_{j,\pm}^k(t,\tau,x,y) = \pm \tau(t - \Psi_j^k(x,y))$ parametrize connected components of $\Gamma_\pm^j$ lying over an open neighborhood of the diagonal of the boundary, we can compute the principal term in the asymptotic expansion of $a_{j,\pm}$ appearing in the expression \eqref{osc int}. Since there are $16$ phases for each $j$, corresponding to $1 \leq k \leq 8$ and $\pm \tau > 0$, \eqref{osc int} should actually be a sum of $16$ oscillatory integrals:
\begin{align}\label{Lagrangian distribution}
S_j(t,x,y) = \sum_{\pm} \sum_{k = 1}^8  \int_0^\infty e^{i \Theta_{j, \pm}^k(t,\tau, x, y)} a_{j,k,\pm}(\tau, x, y) d\tau.
\end{align}
Recalling formula \eqref{principal symbol of an FIO} for the principal symbol of a Fourier integral operator, we must compute the Gelfand-Leray form on the critical set $d_\tau \Theta_{j, \pm}^k = 0$. The Leray measure coming from the Gelfand-Leray form is coordinate invariant and it is ultimately more convenient to first introduce boundary normal coordinates, which we now describe following the presentation in \cite{UhlmannLee1989}.
\\
\\
Fix a point $p \in \d \Omega$. For each $q \in \d \Omega$ near $p$, we denote by $\gamma_q(\mu)$ the unit speed geodesic with initial condition conormal to $\d \Omega$ at $q$ and inward pointing. Now denote by $\phi$ the boundary coordinate which parametrizes $\d \Omega$ with respect to arclength, such that $p$ is given by $\phi = 0$. The coordinate $\phi$ can be extended smoothly inside $\Omega$ so that it is constant along $\gamma_q(\mu)$ for every fixed $q$ near $p$. For $\epsilon > 0$ sufficiently small, $(\mu,\phi)$ is then a smooth coordinate system in an $\epsilon$ neighborhood of $p \in \overline{\Omega}$. In these coordinates, the Euclidean metric is locally given by the warped product
$$
g_{\text{Eucl.}} = d\mu^2 + f(\mu,\phi) d\phi^2,
$$
where $f$ is a locally defined function which is smooth up to the boundary. In particular, $f$ satisfies $f(0,\phi) = 1$. This coordinate system is convenient near the boundary because conformal multiples the vector fields $\d/\d \mu$ and $\d/\d \phi$ extend the orthogonal and and tangential gradients respectively in a tubular neighborhood of the boundary. As the canonical relations $\Gamma_{\pm}^j$ involve both $x$ and $y$ variables, we change $x$ variables to $(\mu, \phi)$ and $y$ variables to $(\nu, \theta)$ according to the procedure described above. In \cite{Vig18}, elliptical polar coordinates were instead used of boundary normal coordinates to acheive a similar decomposition into normal and tangential vector fields.
\\
\\
Without loss of generality, we use boundary normal coordinates instead of Euclidean coordinates in the domain of $\Psi_j^k$ from here on. We now compute the Gelfand-Leray form in boundary normal coordinates.
\begin{lemm}\label{GL Form and Critical set}
The canonical relation of each operator in \eqref{Lagrangian distribution} is parametrized in boundary normal coordinates by
$$
\Gamma_\pm^{j,k} = \{(t,\tau, \mu, \phi, \xi, \nu, \theta, \eta) = (\Psi_j^k, \tau, \mu, \phi, \mp \tau d_{\mu,\phi} \Psi_j^k, \nu, \theta, \pm \tau d_{\nu, \theta} \Psi_j^k) \}.
$$
The Gelfand-Leray form on $C_{\Theta_{j, \pm}^k}$ is given by
\begin{align*}
dC_{\Theta_{j, \pm}^k} =\mp f(\mu, \phi) f(\nu, \theta) d\tau \wedge d\mu \wedge d\phi \wedge d\nu \wedge d\theta.
\end{align*}
\end{lemm}
\begin{proof}
The first claim follows directly from Lemma \ref{parameterization of cannonical relations} and the observation that
$$
d_\tau {\Theta_{j, \pm}^k} = \pm (t - \Psi_j^k) = 0
$$
on the critical set $C_{\Theta_{j, \pm}^k}$. From this, it is clear that $(\tau, \mu, \phi, \nu, \theta)$ form a smooth coordinate system on $C_{\Theta_{j, \pm}^k}$. The Gelfand-Leray form is uniquely defined on $C_{\Theta_{j, \pm}^k}$ by the condition
\begin{align}\label{GLF}
d \left(d_\tau{\Theta_{j, \pm}^k}\right) \wedge dC_{\Theta_{j, \pm}^k} = f(\mu, \phi) f(\nu, \theta) dt \wedge d \tau \wedge d(\mu, \phi) \wedge d(\nu, \theta),
\end{align}
where the righthand side of \eqref{GLF} coincides with the Euclidian volume form on $\mathbb{R}_t \times \mathbb{R}_\tau \times \Omega \times \Omega$. Hence,
$$
dC_{\Theta_{j, \pm}^k} =\mp f(\mu, \phi) f(\nu, \theta) d\tau \wedge d\mu \wedge d\phi \wedge d\nu \wedge d\theta,
$$
as claimed.
\end{proof}
\noindent We now change variables and use Lemma \ref{GL Form and Critical set} to compute the principal symbol of \eqref{Lagrangian distribution}. Dropping the $j,k$ indices on $\Psi_j^k$ in place of differentiation, we have
\begin{align*}
dt &= \Psi_{\mu} d\mu + \Psi_{\phi}d\phi +\Psi_{\nu} d\nu + \Psi_{\theta}d\theta,\\
dy &= f(\nu, \theta) d\nu \wedge d\theta,\\
d{\eta_1} &= \Psi_{\nu} d\tau + \tau( \Psi_{\mu \nu} d\mu + \Psi_{\nu \phi} d\phi + \Psi_{\nu \nu} d\nu + \Psi_{\nu \theta} d\theta),\\
d\eta_2 &= \Psi_{\theta} d\tau +  \tau( \Psi_{ \theta\mu} d\mu + \Psi_{\theta \phi} d\phi + \Psi_{ \theta \nu} d\nu + \Psi_{ \theta \theta} d\theta  )
\end{align*}
on the canonical relation $\Gamma_{\pm}^{j,k}$. Wedging all these terms together, we see that
\begin{align}\label{Aj}
\begin{split}
dt \wedge dy \wedge d\eta = & \tau (f(\nu, \theta)  ( \Psi_{\mu} \Psi_{ \theta} \Psi_{\nu \phi} + \Psi_{\phi} \Psi_{\nu} \Psi_{ \theta \mu}\\
&- \Psi_{\mu}\Psi_{\nu} \Psi_{\theta \phi} - \Psi_{\phi} \Psi_{\theta} \Psi_{\nu \mu} )) d\tau \wedge d\mu \wedge d\phi \wedge d\nu \wedge d\theta.
\end{split}
\end{align}
\begin{def1}
Denote by $A_j^k(\mu,\phi, \nu, \theta)$ the functions
$$
f(\mu, \phi)  ( \Psi_{\mu} \Psi_{ \theta} \Psi_{\nu \phi} + \Psi_{\phi} \Psi_{\nu} \Psi_{ \theta \mu} - \Psi_{\mu}\Psi_{\nu} \Psi_{\theta \phi} - \Psi_{\phi} \Psi_{\theta} \Psi_{\nu \mu} ),
$$
where each $\Psi = \Psi_j^k$ depends implicitly on $j, k, \mu,\phi,\nu$ and $\theta$.
\end{def1}
\noindent On each of the canonical relations $\Gamma_\pm^{j,k}$ (cf. Lemma \ref{parameterization of cannonical relations}), equation \eqref{Aj} implies that
\begin{align}\label{symbol D}
e_\pm = (-1)^j \frac{|dt \wedge dy \wedge d\eta|^{1/2}}{2\tau i} = (-1)^j \frac{\pm 1}{2 |\tau|^{1/2} i}|A_j^k(\mu, \phi, \nu, \theta) dC_{\Theta_{j, \pm}^k}|^{1/2}
\end{align}
for Dirichlet boundary conditions and
\begin{align}\label{symbol N}
e_\pm = \frac{|dt \wedge dy \wedge d\eta|^{1/2}}{2\tau i} = \frac{\pm 1}{2 |\tau|^{1/2} i}|A_j^k(\mu, \phi, \nu , \theta) dC_{\Theta_{j, \pm}^k}|^{1/2}
\end{align}
for Neumann boundary conditions. As a result, we conclude the following description of the operators $S_j$.
\begin{theo}\label{parametrix with sum}
Microlocally near $\Gamma_\pm^j$, the following oscillatory integral is a parametrix for $S_j(t,x,y)$ with Dirichlet boundary conditions in an open neighborhood of $\Delta \d \Omega \subset \Omega \times \Omega$:
\begin{align*}
S_j(t)  = (-1)^j \sum_{\pm} \sum_{k = 1}^8 \int_{0}^\infty e^{\pm i\tau(t - \Psi_j^k)} \frac{\pm 1 }{2|\tau|^{1/2} i} |A_j^k|^{1/2} d\tau + L.O.T.
\end{align*}
Here, L.O.T. denotes lower order Lagrangian distributions, using the convention \eqref{order of a lagrangian distribution}.
\end{theo}
\begin{def1}
We define the operators $S_{j,\pm}^k(t)$ appearing in Theorem \ref{parametrix with sum} by
$$
S_{j,\pm}^k(t) =  (-1)^j \sum_{\pm} \sum_{k = 1}^8 \int_{0}^\infty e^{\pm i\tau(t - \Psi_j^k)} \frac{\pm 1 }{2|\tau|^{1/2} i} |A_j^k|^{1/2} d\tau.
$$
\end{def1}
\noindent Together with Theorem \ref{8 orbit lemma} and the choice of $j_0$, $U_j$ in Remark \ref{Uj}, we can immediately derive Theorem \ref{HRP} from Theorem \ref{parametrix with sum}. To see this, note that
$$
E(t) = \frac{d}{dt} S(t).
$$
On a symbolic level, this implies that the symbol of $E(t)$ is $\tau e_\pm$, with $e_\pm$ given by Proposition \ref{symbol prop}. This explains the order of $b_{j,k,\pm} \in S^{1/2}( U_j \times \mathbb{R}^1)$ appearing in Theorem \ref{HRP}. Each $(\pm)$ branch of the propagators $S_{j,\pm}^k$ should be multiplied by Maslov factors $e^{\pm i \pi \sigma_j/4}$. We will compute $\sigma_j = 1$ in Section \ref{Maslov factors}, which completes the proof of Theorem \ref{HRP}. In the next section, we will use this explicit oscillatory integral representation to compute the wave trace.

\section{Computing the wave trace}\label{Computing the wave trace}
\noindent In this section, we use the parametrix in Theorem \ref{HRP}, or equivalently Theorem \ref{parametrix with sum}, to compute an integral formula for the wave trace. Formally, the wave trace is the Fourier transform of the spectral measure $\sum_j \delta(\lambda - \lambda_j)$, where $\{\lambda_j^2\}$ are the Dirichlet (resp. Neumann or Robin) eigenvalues of $-\Delta$ on $\Omega$. This is a distribution of the form
\begin{align}\label{eigenvalue wave trace}
\sum_j e^{i t \lambda_j},
\end{align}
which can be seen to be weakly convergent by Weyl's law on the asymptotic distribution of Laplace eigenvalues (see \cite{Irvii16}, \cite{Zayed2004}). The connection between this distribution and the wave equation lies in the fact that \eqref{eigenvalue wave trace} is actually the trace of $e^{i t \sqrt{- \Delta}}$, the propagator associated to the half wave operator $(\d_t - i \sqrt{-\Delta})$. Since such a unitary operator is not trace class, we mean that for any Schwartz function $\phi$, the regularized operators
$$
\int_{-\infty}^\infty \phi(t) e^{it \sqrt{-\Delta}}dt
$$
are of trace class and have trace
$$
\sum_j \int_{-\infty}^\infty \phi(t) e^{it \lambda_j}dt.
$$
The same holds for the even and odd wave operators
\begin{align}
E(t) = \cos t\sqrt{- \Delta}, \qquad S(t) = \frac{\sin t \sqrt{-\Delta}}{\sqrt{-\Delta}},
\end{align}
and we consider these as they appear more naturally in Chazarain's parametrix (cf. Section \ref{Chazarain's Parametrix}) and $E(t)$ solves the initial boundary value problem \eqref{wave equation}.
\\
\\
Recall that $\chi_1(t)$ is a cutoff function near the lengths $[t_j,T_j]$ of geodesic loops and $\Omega$ is assumed to satisfy the noncoincidence condition \eqref{NCC}. Modulo Maslov factors and a smooth error term, we have
\begin{align}\label{localized wave trace}
\chi_1(t) \text{Tr} E &= \int_{\Omega} E_j(t,x,x) dx,
\end{align}
where $E_j = \frac{\d}{\d t} S_j$ and $S_j$ is the $j$ bounce wave appearing in Chazarain's parametrix \eqref{Chazarain sum}.
\\
\\
As our parametrix in Theorem \ref{parametrix with sum} is only valid near the boundary, we want to write the wave trace as an integral over the boundary. There are two ways to do this. The first involves writing $E_j = \Delta B_j$ for a specially chosen distribution $B_j$ so that formally, we may integrate by parts to obtain an integral of $\nabla B_j$ in terms of $E_j$ over the boundary. The disadvantage of this approach is nonlocality of $\Delta^{-1}$, which makes it difficult to compute $B_j$ explicitly. The second way is to use Hadamard type variational formulas for the wave trace, which were derived in \cite{HeZe12} and \cite{Vig18}. As shown in those papers, taking the variation of the wave trace automatically includes an integration by parts which puts the integral \eqref{localized wave trace} on the boundary. The advantage of this approach is that the formulas are quite explicit.

\subsection{Rescaling techniques}
Our idea is to rescale the domain $\Omega$ by dilating about an interior point. Setting $\Omega_r = r \Omega$ for $0 \leq r \leq 1$, the domain shrinks to a point when $r = 0$. Denote by $\Delta_r$ the Laplacian on $\Omega_r$ with Dirichlet boundary conditions. Formally, one then has
\begin{align}\label{intuition hadamard}
\int_{0^+}^1 \frac{\d}{\d r} \text{Tr} \cos t \sqrt{-\Delta_{r}} dr = \text{Tr} \cos t\sqrt{-\Delta_{1}} - \lim_{\epsilon \to 0} \text{Tr} \cos t \sqrt{-\Delta_{\epsilon}}.
\end{align}
Here, the integral bounds $(0^+,1]$ on the lefthand side of \eqref{intuition hadamard} correspond to a distributional limit, where the domain of integration is taken to be $(\epsilon, 1]$ as $\epsilon \to 0^+$. Equation \eqref{intuition hadamard} is essentially the fundamental theorem of calculus. The quantity $\d_r \cos t\sqrt{-\Delta_{\Omega_r}}$ is given explicitly by variational formulas for the wave trace in \cite{Vig18} and \cite{HeZe12}. As $r \to 0^+$, the eigenvalues of $-\Delta_r$ increase rapidly and monotonically to $+ \infty$. This contributes highly oscillatory terms to the wave trace \eqref{eigenvalue wave trace}, which tend weakly to zero with $r$. We will in fact only scale down to $r = r_0$ rather than $r = 0$, for $r_0$ sufficiently small at the expense of smooth errors. This is explained in detail in Section \ref{Integrating the principle symbol}. In this section, we derive more explicit integral formulas for the lefthand side of equation \eqref{intuition hadamard}. In Section \ref{Principle symbol computation}, we will then justify these computations in more detail and show that they indeed give the trace of $\cos t \sqrt{-\Delta_1}$.
\\
\\
As $\Delta$ is translation invariant, we may assume without loss of generality that the origin $O$ is an interior point of $\Omega$. For each $r \in (0,1]$, consider the dilated domain $\Omega_r = r \Omega$ with respect to $O$. Define wave propagators on the resclaed domains $\Omega_r$ by
$$
E_r(t) = \cos t \sqrt{- \Delta_r} \qquad \qquad S_r(t) = \frac{\sin t \sqrt{-\Delta_r}}{\sqrt{-\Delta_r}}.
$$
The formulas in \cite{HeZe12} and \cite{Vig18} show that
\begin{align}\label{Hadamard reference}
\d_r \text{Tr} \cos t \sqrt{-\Delta_r} = \begin{cases}
\frac{t}{2} \int_{\d \Omega_r} \nabla_1^\perp \nabla_2^\perp S(t,q,q) \dot{\rho}\,dq & \text{Dirichlet BC}\\
\frac{-t}{2} \int_{\d \Omega_r} (\nabla_1^T \nabla_2^T  + \Delta_2) S(t,q,q) \dot{\rho}\,dq & \text{Neumann BC},
\end{cases}
\end{align}
where $\nabla_{1,2}^\perp$ are the unit orthogonal gradients on $\d \Omega_r$ and $\dot{\rho}$ is the normal component of the variational vector field $\d_r$, which will be computed in Section \ref{Principle symbol computation}. Note that the integrals in \eqref{Hadamard reference} are over the diagonal of the boundary. We now discuss how to calculate these integrals more precisely.

\subsection{Wave propagators near the boundary}
As we are localizing the wave trace near lengths of geodesic loops with $j$ reflections, it turns out we only need to insert a select few of the operators from Theorem \ref{parametrix with sum} into the formula \eqref{Hadamard reference}.
\begin{lemm}\label{Sj lemma}
Modulo Maslov factors and lower order distributions, the radially differentiated even wave trace localized near $[t_j,T_j]$ is given by
\begin{align*}
\sum_\pm \int_{\d \Omega} \frac{t}{2} \nabla_1^\perp \nabla_2^\perp (&S_{j-1,\pm}^1  +S_{j-1,\pm}^5 + S_{j,\pm}^2 + S_{j,\pm}^3 \\
& + S_{j,\pm}^6 + S_{j,\pm}^7 + S_{j + 1,\pm}^4 +  S_{j +1,\pm}^8 )(t,q,q) \dot{\rho}\,dq
\end{align*}
in the Dirichlet case and
\begin{align*}
\sum_\pm \int_{\d \Omega} \frac{t}{2} (-\nabla_1^T \nabla_2^T - \Delta_2)(&S_{j-1,\pm}^1  +S_{j-1,\pm}^5 + S_{j,\pm}^2 + S_{j,\pm}^3 \\
& + S_{j,\pm}^6 + S_{j,\pm}^7 + S_{j + 1,\pm}^4 +  S_{j +1,\pm}^8 )(t,q,q) \dot{\rho}\,dq
\end{align*}
in the Neumann case.
\end{lemm}
\noindent For completeness, we repeat the proof derived in \cite{Vig18} as it contains substantial geometric insight.
\begin{proof}
For the localized wave trace, we only need to consider orbits which contribute to the singularities in $[t_j,T_j]$. Recall that for positive time, Theorem \ref{8 orbit lemma} gives $8$ orbits connecting $x$ to $y$ in $j$ reflections and approximately one rotation. These orbits coalesce into geodesic loops as $(x, y) \to \Delta \d \Omega$. However, as the orbits coalesce within various configurations, not all of the limiting orbits will have $j$ reflections. As $\Omega$ satisfies the noncoincidence condition \eqref{NCC}, only the limiting geodesic loops of rotation number $1/j$ will contribute to the wave trace near $[t_j,T_j]$. See Figure \ref{configurations} for visualizing the geometric arguments which follow. As $(x, y) \to \Delta \d \Omega$, the two orbits of $j$ reflections in $TT$ configuration $(k = 1,5)$ converge geometrically to a loop of $j+1$ reflections. The additional vertex appears at the boundary point where $x$ and $y$ coalesce. Similarly, the two $NN$ orbits of $j$ reflections $(k = 4,8)$ can be seen to converge to loops of $j-1$ reflections. In this case, the first and last moments of reflection at the boundary converge to a single impact point. The four orbits of $j$ reflections in $TN$ $(k = 2,6)$ and $NT$ $(k = 3,7)$ configurations preserve exactly $j$ reflections in the limit. Hence, when $x, y \in \text{int} \Omega$ converge to $\Delta \d \Omega$, only $4$ of the $8$ orbits contribute to geodesic loops of $j$ reflections. However, in the same limit, two additional $TT$ orbits of $j-1$ reflections converge to a loop of $j - 1 + 1 = j$ reflections. Similarly, two $NN$ orbits of $j+1$ reflections converge to a loop of $j+1 -1 = j$ reflections. Any other orbit from $x$ to $y$ with strictly less than $j-1$ or strictly more than $j+1$ reflections at the boundary cannot converge to a loop of $j$ reflections. As we have localized the wave trace near the isolated set of lengths $[t_j,T_j]$ and $\Omega$ satisfies the noncoincidence condition \eqref{NCC}, only the $4 + 2 + 2 = 8$ orbits which converge geometrically to a loop of exactly $j$ reflections will contribute to singularities here. All additional orbits contribute smooth errors to the wave trace in a small neighborhood of $[t_j,T_j]$.
\\
\\
It should also be clarified that although the parametrices $S_j(t,x,y)$ are constructed in the interior, we can in fact extend them continuously to the diagonal of the boundary and this extension coincides with that of the true propagator $S(t,x,x)$ ($x \in \d \Omega$) modulo lower order terms. Both propagators agree up to lower order Lagrangian distributions in the interior, microlocally near the canonical relations $\Gamma_\pm^j$. The explicit oscillatory integral representation for each $S_j(t,x,y)$ in fact shows that they extend continuously up to the boundary and its diagonal, since the functions $\Psi_j^k(x,y)$ do. The true wave kernel $S(t,x,y)$ also extends continuously up to the boundary as a family of distributions. To see this, note that
$$
S(t,x,y) = \sum_j \frac{\sin ( t \lambda_j)}{\lambda_j} \psi_j(x) \overline{\psi_j(y)},
$$
where $(\psi_j)_{j =1}^\infty$ is an $L^2$ orthonormal basis of Dirichlet or Neumann eignenfunctions corresponding to eigenvalues $(\lambda_j^2)_{j = 1}^\infty$. Multiplying by a test function $\phi(t)$ and integrating by parts $4k$ times, we see that
\begin{align}\label{wave kernel}
\int_{-\infty}^\infty S(t,x,y) \phi(t) dt = \int_{-\infty}^\infty \sum_j \frac{\sin (t \lambda_j)}{\lambda_j^{4k+1}} \psi_j(x) \overline{\psi_j(y)} \d_t^{4k}  \phi(t) dt.
\end{align}
Combining Weyl's law on the asymptotic growth of $\lambda_j$ (see \cite{Zayed2004}, \cite{Irvii16}) and H\"ormander's $L^\infty$ bounds for eigenfunctions (see \cite{HormanderLinfinity}), we see that the integrand in \eqref{wave kernel} can be made absolutely convergent for $k$ sufficiently large. An application of the dominated convergence theorem then shows that \eqref{wave kernel} is actually smooth in $x,y$, so $S(t,x,y)$ has a smooth extension to the diagonal of the boundary as a distribution in $t$. In particular, both distributions agree up to lower order terms microlocally near the fibers of $\Gamma_\pm^j$ lying over diagonal of the boundary, as was required for the variational trace formula.
\end{proof}

\begin{def1}\label{psi j coincide on boundary}
As shown in the proof of Lemma \ref{Sj lemma} above, for each $j$, there exist $8$ limiting trajectories which converge geometrically to geodesic loops of exactly $j$ reflections. We denote the set of these trajectories by $\mathcal{G}_j(x,y)$ and say that $\gamma_{m,k} \in \mathcal{G}_j$ if $\gamma_{m,k}$ makes $m = j-1, j$ or $j+1$ reflections at the boundary and corresponds to the length functional $\Psi_m^k$. The length functionals $\Psi_j^2,\Psi_j^3 , \Psi_j^6 , \Psi_j^7 , \Psi_{j + 1}^4 , \Psi_{j +1}^8 , \Psi_{j-1}^1$ and $\Psi_{j-1}^5$ corresponding to orbits in $\mathcal{G}_j$ coincide for $x =y \in \d \Omega$ on the diagonal. We denote their common value by $\Psi_j(x,x)$, which is the $j$-loop function..
\end{def1}
\noindent As we obtained a rather explicit formula for $S_j(t,x,y)$ in Theorem \ref{parametrix with sum}, it now remains to differentiate the kernels $S_m^k(t,x,y)$ and substitue them into Lemma \ref{Sj lemma}. Using our oscillatory integral representation for $S_m^k(t,x,y)$ in Theorem \ref{parametrix with sum}, we find that microlocally near $\Gamma_\pm^{m,k}$ and $t \in [t_j,T_j]$, modulo lower order terms in an open neighborhood of the diagonal of the boundary, we have
\begin{align}\label{parametrix with 8 Dirichlet}
\nabla_1^\perp \nabla_2^\perp S_{m,\pm}^k = (-1)^{m+1} \sum_\pm \int_{0}^\infty e^{\pm i\tau(t - \Psi_m^k)} \frac{\pm (\nabla_1^\perp \Psi_m^k)(\nabla_2^\perp \Psi_m^k) |\tau|^{3/2}}{2 i} |A_m^k|^{1/2} d\tau
\end{align}
for Dirichlet boundary conditions.
The formula for Neumann boundary conditions is less succinct and from here on, we only study Dirichlet boundary conditions. We have only written the terms coming from $\nabla_{1,2}^\perp$ falling on the exponential in equation \eqref{parametrix with 8 Dirichlet}, as all other terms don't contribute positive powers of $\tau$ and can be regarded as lower order in the singularity expansion. The operators $\nabla_{1,2}^\perp$ in the integrand of \eqref{parametrix with 8 Dirichlet} are conformal multiples of the vector fields $\frac{\d}{\d \mu}$ and $\frac{\d}{\d \nu}$ coming from boundary normal coordinates.
\\
\\
As equation \eqref{Hadamard reference} tells us that the variation of the wave trace is given by integrating the differentiated sine kernels over the diagonal of the boundary, we want to understand the restriction of \eqref{parametrix with 8 Dirichlet} to the boundary. We already noted that the $j$-loop function is well defined for $j$ sufficiently large.
The differentiated kernels in equation \eqref{parametrix with 8 Dirichlet} also have factors of $A_j^k$ and $\nabla_{1,2}^\perp \Psi_j^k$ in the integrand. We now discuss how to extend these derivatives of $\Psi_j^k$ to the diagonal of the boundary in a similar manner. In the proof of Lemma \ref{parameterization of cannonical relations} (cf. \cite{Vig18}), the $x$ and $y$ gradients of the functions $\Psi_j^k$ are shown to be:
\begin{align}\label{x y gradients}
d_x \Psi_j^k = \frac{x - q_1^k}{|x - q_1^k|}, \quad d_y \Psi_j^k = \frac{y - q_j^k}{|y - q_j^k|}.
\end{align}
Geometrically, these are the incident and reflected outgoing unit directions of the corresponding billiard trajectories at $x$ and $y$. The expression $\nabla_1^\perp \Psi_j^k$ in \eqref{parametrix with 8 Dirichlet} can easily be seen to be $\pm \sin \omega_{j,1}^k$, where $\omega_{j,1}^k$ is the angle made between the initial link of the billiard trajectory and the oriented tangent line to the level set of the distance function on which $x$ lies (as in the folliation used in the proof of Theorem \ref{8 orbit lemma}). We use the positively oriented tangent line for $1 \leq k \leq 4$ and the negatively oriented tangent line for $5 \leq k \leq 8$. The sign $\pm$ depends on the $TT, TN, NT$ or $NN$ configuration of the corresponding orbit. Similarly, $\nabla_2^\perp \Psi_j^k = \pm \sin \omega_{j,2}^k$, where $\omega_{j,2}^k$ is the angle made between the final link of the billiard trajectory and the positively $(1 \leq k \leq 4)$ or negatively $(5 \leq k \leq 8)$ oriented tangent line to the distance curve on which $y$ lies. As $x,y \to  \Delta \d \Omega$, the absolute value of the angles associated to trajectories in the $\mathcal{G}_j$ converge to the initial and final angles of reflection of the unique limiting geodesic loop. We are careful to point out that only the absolute values of the angles converge, since for example, the final angles of reflection at $y$ associated to orbits in $TN$ and $NN$ configurations actually converge to the \textit{negative} of the final angle in the limiting trajectory. This phenomenon will be important for avoiding cancellation in the equation \eqref{D1} of Section \ref{Integrating the principle symbol}, where we add together all wave kernels and integrate them over the diagonal to obtain a trace formula. All limiting loops are automatically in $TT$ configuration.
\begin{lemm}\label{Aj w derivatives}
On the diagonal of the boundary, the factors $A_{j-1}^1, A_{j-1}^5 A_j^2, A_j^3, A_j^6, A_j^7, A_{j+1}^4$ and $A_{j+1}^8$ corresponding to orbits in $\mathcal{G}_j$ coincide up to a sign. We denote their common (absolute) value by $|A_j|$, which at a point $(0,\phi, 0, \phi)$ in boundary normal coordinates satisfies
\begin{align*}
|A_j(0, \phi, 0, \phi)| =\frac{1}{\sin \omega_{j,2}}  \left| \frac{\d \omega_{j,1}}{\d \theta} \right|.
\end{align*}
Here, $\omega_{j,1}$ and $\omega_{j,2}$ are the initial and final angles of incidence respectively made by the unique counterclockwise parametrized geodesic loop of rotation number $1/j$ based at $(0, \phi, 0, \phi) \in \Delta \d \Omega$.
\end{lemm}
\begin{proof}
Let $(m,k)$ denote an admissable pair of indices corresponding to an orbit $\gamma_{m,k} \in \mathcal{G}_j$. Recall the notation in the proof of Lemma \ref{parameterization of cannonical relations} (cf. \cite{Vig18}), where we described a billiard trajectory by the point $(x,q,y) \in \Omega \times \d \Omega^m \times \Omega$. Let us first assume $1 \leq k \leq 4$. If $x \in \Omega$, we denote the positive angle between $q_1 - x$ and the positively oriented tangent line to the leaf of the folliation by distance curves on which $x$ lies by $\omega_{j,1}^k$ (cf. Section \ref{Proof of 8 orbit lemma}). Similarly, if $y \in \Omega$, let us also denote the positive angle between $y- q_j$ and the positively oriented tangent line to the distance curve on which $y$ lies by $\omega_{j,2}^k$. Note that $\omega_{j,1}^k$ and $\omega_{j,2}^k$ depend on $x, y, j$ and $k$. Using the warped product structure of boundary normal coordinates $(\mu, \phi, \nu, \theta)$, we have
\begin{align}\label{BNCG}
\begin{split}
\nabla_x^T &= {(f)}^{-1/2} \frac{\d}{\d\phi}, \qquad
\nabla_y^T = {(f)}^{-1/2} \frac{\d}{\d\theta},\\
\nabla_x^\perp &=  \frac{\d}{\d\mu}, \qquad \qquad \,\,\,\,\,\,\,\,
\nabla_y^\perp =  \frac{\d}{\d\nu}.
\end{split}
\end{align}
Equation \eqref{x y gradients}, which was derived from the proof of Lemma \ref{parameterization of cannonical relations} in \cite{Vig18}, then shows that
\begin{align}\label{perp grads}
\begin{split}
\nabla_x^\perp \Psi_j^k(x,y) = \begin{cases}
- \sin \omega_{j,1}^k, & k = 1,2\\
\sin \omega_{j,1}^k, & k = 3,4,
\end{cases} \qquad
\nabla_y^\perp \Psi_j^k(x,y) = \begin{cases}
- \sin \omega_{j,2}^k, & k = 1,3\\
\sin \omega_{j,2}^k, & k = 2,4,
\end{cases}
\end{split}
\end{align}
and
\begin{align}\label{tangential grads}
\nabla_x^T(x,y) = -\cos \omega_{j,1}^k, \quad (1 \leq k \leq 4) \qquad \qquad \nabla_y^T = \cos \omega_{j,2}^k, \quad (1 \leq k \leq 4).
\end{align}
Combining \eqref{BNCG}, \eqref{perp grads}, and \eqref{tangential grads}, we obtain
\begin{align}\label{first derivatives}
\begin{split}
\nabla_{(\mu,\phi)} \Psi_j^k(\mu, \phi, \nu, \theta) &= \begin{cases}
(-\sin \omega_{j,1}^k, -f^{1/2}(\mu, \phi) \cos \omega_{j,1}^k ), & k = 1\\
(-\sin \omega_{j,1}^k, -f^{1/2}(\mu, \phi) \cos \omega_{j,1}^k ), & k = 2\\
(\sin \omega_{j,1}^k, -f^{1/2}(\mu, \phi) \cos \omega_{j,1}^k ), & k = 3\\
(\sin \omega_{j,1}^k, -f^{1/2}(\mu, \phi) \cos \omega_{j,1}^k ), & k = 4
\end{cases}\\
\nabla_{(\nu,\theta)} \Psi_j^k(\mu, \phi, \nu, \theta) &= \begin{cases}
(-\sin \omega_{j,2}^k, f^{1/2}(\nu, \theta) \cos \omega_{j,2}^k ),& k = 1\\
(\sin \omega_{j,2}^k, f^{1/2}(\nu, \theta) \cos \omega_{j,2}^k ), & k = 2\\
(-\sin \omega_{j,2}^k, f^{1/2}(\nu, \theta) \cos \omega_{j,2}^k ),& k = 3\\
(\sin \omega_{j,2}^k, f^{1/2}(\nu, \theta) \cos \omega_{j,2}^k ), & k = 4.
\end{cases}
\end{split}
\end{align}
Using \eqref{first derivatives} to calculate the $(\mu,\phi,\nu,\theta)$ Hessian of $\Psi_j^k$, we find
\begin{align}\label{second derivatives}
\begin{split}
\frac{\d^2 \Psi_j^k}{\d \theta \d \phi}(\mu, \phi,\nu, \theta) &=  f^{1/2}(\mu,\phi) \sin \omega_{j,1}^k \frac{\d \omega_{j,1}^k }{\d \theta}, \qquad (1 \leq k \leq 4)\\
\frac{\d^2 \Psi_j^k}{\d \nu \d \phi}(\mu, \phi,\nu, \theta) &= f^{1/2}(\mu,\phi) \sin \omega_{j,1}^k \frac{\d \omega_{j,1}^k }{\d \nu}, \qquad (1 \leq k \leq 4)\\
\frac{\d^2 \Psi_j^k}{\d \theta \d \mu}(\mu, \phi,\nu, \theta) &= \begin{cases}
- \cos \omega_{j,1}^k \frac{\d \omega_{j,1}^k }{\d \theta}, & k = 1,2\\
\cos \omega_{j,1}^k \frac{\d \omega_{j,1}^k }{\d \theta}, & k = 3,4
\end{cases} \\
\frac{\d^2 \Psi_j^k}{\d \nu \d \mu}(\mu, \phi,\nu, \theta) &= \begin{cases}
- \cos \omega_{j,1}^k \frac{\d \omega_{j,1}^k }{\d \nu}, & k = 1,2\\
\cos \omega_{j,1}^k \frac{\d \omega_{j,1}^k }{\d \nu}, & k = 3,4.
\end{cases}
\end{split}
\end{align}
Inserting \eqref{second derivatives} into the expression \eqref{Aj} for $A_j^k$ in all possible configurations $(1 \leq k \leq 4)$, we find that on the boundary,
\begin{align}\label{new Aj w derivatives}
A_m^k(0, \phi, 0, \theta) =
\begin{cases}
-\cos \omega_{j,2}^k \frac{\d\omega_{j,1}^k}{\d \nu} - \sin \omega_{j,2}^k \frac{\d \omega_{j,1}^k}{\d \theta}, & k = 1\\
- \cos \omega_{j,2}^k \frac{\d\omega_{j,1}^k}{\d \nu} + \sin \omega_{j,2}^k \frac{\d \omega_{j,1}^k}{\d \theta},  & k = 2\\
+ \cos \omega_{j,2}^k \frac{\d \omega_{j,1}^k}{\d \nu} + \sin \omega_{j,2}^k \frac{\d \omega_{j,1}^k}{\d \theta}, & k = 3\\
+ \cos \omega_{j,2}^k \frac{\d\omega_{j,1}^k}{\d \nu} - \sin \omega_{j,2}^k \frac{\d \omega_{j,1}^k}{\d \theta}, & k = 4.
\end{cases}
\end{align}
Before evaluating this expression on the diagonal of the boundary, we differentiate $\omega_{j,1}^k$ in the direction $L = (y-q)/|y-q|$ of the last link to see that
\begin{align*}
0 &= \nabla_L \omega_{j,1}^k\\
&= \begin{cases}
\cos \omega_{j,2}^k \nabla_y^T \omega_{j,1}^k - \sin \omega_{j,2}^k \nabla_y^\perp \omega_{j,1}^k, & k = 1,3\\
\cos \omega_{j,2}^k \nabla_y^T \omega_{j,1}^k + \sin \omega_{j,2}^k \nabla_y^\perp \omega_{j,1}^k, & k = 2,4
\end{cases}\\
&= \begin{cases}
\frac{1}{f^{1/2}} \cos \omega_{j,2}^k  \frac{\d \omega_{j,1}^k}{\d \theta} - \sin \omega_{j,2}^k  \frac{\d \omega_{j,1}^k}{\d \nu}, & k = 1,3\\
\frac{1}{f^{1/2}} \cos \omega_{j,2}^k  \frac{\d \omega_{j,1}^k}{\d \theta} + \sin \omega_{j,2}^k  \frac{\d \omega_{j,1}^k}{\d \nu}, & k = 2,4.
\end{cases}
\end{align*}
This implies that
\begin{align}\label{d omega d nu}
\sqrt{f(\nu, \theta)} \frac{\d \omega_{j,1}^k}{\d \nu} = \begin{cases}
+\cot \omega_{j,2}^k \frac{\d \omega_{j,1}^k}{\d \theta}, & k =1,3\\
- \cot \omega_{j,2}^k \frac{\d \omega_{j,1}^k}{\d \theta}, & k = 2,4.
\end{cases}
\end{align}
Inserting the formula \eqref{d omega d nu} into \eqref{new Aj w derivatives} and evaluating on $\d \Omega$, we find that
\begin{align*}
|A_m^k(0, \phi, 0, \theta)| =  \left(\frac{\cos^2 \omega_{j,2}^k}{\sin \omega_{j,2}^k} + \frac{\sin^2 \omega_{j,2}^k}{\sin \omega_{j,2}^k}  \right) \left| \frac{\d \omega_{j,1}^k}{\d \theta} \right|  = \frac{1}{\sin \omega_{j,2}^k} \left| \frac{\d \omega_{j,1}^k}{\d \theta} \right|.
\end{align*}
Note that for $1 \leq k \leq 4$ and $x, y$ on the boundary, both $\omega_{j,1}^k$ and $\omega_{j,2}^k$ are independent of $k$. They coincide with the inital and final angles of the unique orbit connecting $x$ to $y$ in $j$ reflections and approximately one counterclockwise rotation. As the orbits corresponding to $5 \leq k \leq 8$ can be viewed as counterclockwise orbits in the reflected domain, we instead defined $\omega_{j,1}^k$ and $\omega_{j,2}^k$ to be the angles made between the initial and final links and the \textit{negatively} oriented tangent lines to the distance curves on which $x$ and $y$ respectively lie. With this convention, we see by symmetry that the roles of $\omega_{j,1}^k$ and $\omega_{j,2}^k$ are interchanged:
\begin{align}\label{first derivatives k58}
\begin{split}
\nabla_{(\mu,\phi)} \Psi_j^k(\mu, \phi, \nu, \theta) &= \begin{cases}
(-\sin \omega_{j,1}^k, f^{1/2}(\mu, \phi) \cos \omega_{j,1}^k ), & k = 5\\
(-\sin \omega_{j,1}^k, f^{1/2}(\mu, \phi) \cos \omega_{j,1}^k ), & k = 6\\
(\sin \omega_{j,1}^k, f^{1/2}(\mu, \phi) \cos \omega_{j,1}^k ), & k = 7\\
(\sin \omega_{j,1}^k, f^{1/2}(\mu, \phi) \cos \omega_{j,1}^k ), & k = 8
\end{cases}\\
\nabla_{(\nu,\theta)} \Psi_j^k(\mu, \phi, \nu, \theta) &= \begin{cases}
(-\sin \omega_{j,2}^k, -f^{1/2}(\nu, \theta) \cos \omega_{j,2}^k ),& k = 5\\
(\sin \omega_{j,2}^k, -f^{1/2}(\nu, \theta) \cos \omega_{j,2}^k ), & k = 6\\
(-\sin \omega_{j,2}^k, -f^{1/2}(\nu, \theta) \cos \omega_{j,2}^k ),& k = 7\\
(\sin \omega_{j,2}^k, -f^{1/2}(\nu, \theta) \cos \omega_{j,2}^k ), & k = 8.
\end{cases}
\end{split}
\end{align}
Replacing \eqref{first derivatives} by \eqref{first derivatives k58}, parallel computations to those above then show that for $5 \leq k \leq 8$,
\begin{align*}
|A_m^k(0,\phi,0,\phi)| = \frac{1}{\sin \omega_{j,2}^k} \Big|\frac{\d \omega_{j,1}}{\d \theta}\Big|.
\end{align*}
Note that for $5 \leq k \leq 8$, $\omega_{j,1}^{k-4} - \omega_{j,2}^k \to 0$ as $x, y \to \d \Omega$. Similarly, $\omega_{j,2}^{k-4} - \omega_{j,1}^k \to 0$ as $x, y\to \d \Omega$. Hence,
\begin{align}\label{k58}
|A_m^k(0,\phi,0,\phi)| = \frac{1}{\sin \omega_{j,1}} \Big| \frac{\d \omega_{j,2} }{ \d \phi} \Big|,
\end{align}
for $5 \leq k \leq 8$. Note that the righthand side of \eqref{k58} involves angles at boundary points and does not depend on $k$. Also observe that
\begin{align}\label{Clairaut}
\begin{split}
\frac{\d^2 \Psi_j^k}{\d \phi \d \theta}  &= \pm \sin \omega_{j,2}^k \frac{\d \omega_{j,2}^k}{\d \phi}\\
\frac{\d^2 \Psi_j^k}{ \d \theta \d \phi} &= \pm \sin \omega_{j,1}^k \frac{\d \omega_{j,1}^k}{\d \theta}.
\end{split}
\end{align}
Setting the two equations in \eqref{Clairaut} equal implies that
\begin{align*}
\frac{1}{\sin \omega_{j,1}} \Big| \frac{\d \omega_{j,2} }{ \d \phi} \Big| = \frac{1}{\sin \omega_{j,2}} \Big| \frac{\d \omega_{j,2} }{ \d \phi} \Big|,
\end{align*}
which combined with \eqref{k58}, concludes the proof of the lemma.
\end{proof}
\begin{rema}
In \cite{Vig18}, complete integrability of the ellipse actually implies that $\omega_{j,1} = \omega_{j,2}$ for elliptical billiards i.e. every geodesic loop is in fact a periodic orbit. In that case, the angular derivative $\d \omega / \d \theta$ appearing in Lemma \ref{Aj w derivatives} was explicitly calculated, using action angle coordinates and Jacobi elliptic functions. In general, it may be difficult to compute $\d \omega_1/\d \theta$ explicitly despite its relatively simple geometric interpretation.
\end{rema}

\section{Principal symbol computation}\label{Principle symbol computation}
\noindent In this section, we refine the integral formulas obtained in Section \ref{Computing the wave trace} and use them to compute the principal symbol of the wave trace. We also compute Maslov factors on each branch of the wave propagator and provide estimates on the subprincipal terms of the parametrix described in Theorem \ref{parametrix with sum}. These estimates ensure that radially integrating subprincipal terms in the variational wave trace does not contribute to the leading order asymptotics in equation \eqref{intuition hadamard}. Combining these estimates with formulas in Section \ref{Computing the wave trace} then allows us to complete the proofs of Theorem \ref{Main theorem}, Theorem \ref{Rational caustic theorem} and Corollary \ref{Ellipse corollary}.
\\
\\
We begin by evaluating the quantity $\dot{\rho}$, coming from the variational formula \eqref{Hadamard reference}. Let $X(s) = (x(s), y(s))$ be an arclength parametrization of $\d \Omega$. Then, $rX(s) = (rx(s), r y(s))$ parametrizes $\d \Omega_r$ and coincides with a multiple of the variational vector field $\d / \d r$. The unit normal to $\d \Omega_r$ at $(rx,ry)$ is given by $N(s) = (-y'(s), x'(s))$, which is independent of $r$. Hence, $\dot{\rho}(s) = X\cdot N$. Plugging this into \eqref{Hadamard reference} together with \eqref{parametrix with 8 Dirichlet}, we see that locally near $[t_j,T_j]$, modulo Maslov factors and lower order terms, $\d_r \text{Tr} \cos t\sqrt{-\Delta_{\Omega_r}}$ is given by
\begin{align*}
\sum_{\gamma_{m,k} \in \mathcal{G}_j, \pm}\frac{ \pm t  (-1)^{m+1}}{4i} \int_{\d \Omega_r} \int_{0}^\infty \Delta^*( e^{\pm i \tau(t - \Psi_m^k)} |\tau|^{3/2} (\nabla_1^\perp \Psi_m^k)(\nabla_2^\perp \Psi_m^k)|A_m^k|^{1/2}) X\cdot N d\tau ds_r
\end{align*}
(for Dirichlet boundary conditions).
Recall that here, $\Delta^*$ is pullback under the diagonal embedding $\Delta: \d \Omega \to \d \Omega \times \d \Omega$ so that all terms depending on both $x$ and $y$ in integrand above are evaluated at the diagonal of the boundary $\{x = y = s \in \d \Omega\}$. Note that the functions $\Psi_m^k$, $A_m^k$ as well as the operators $\nabla_i^\perp$ and the measure $ds_r$ all depend on $r$. It is easily verified that $\Psi_m^k$ and $ds_r$ are homogeneous of degree one in $r$, $\nabla_{1,2}^\perp \Psi_m^k$ is homogeneous of degree $0$, and $A_m^k$ is homogeneous of degree $-1$. In particular, the orthogonal gradients restricted to the diagonal of the boundary are still given by
\begin{align*}
\nabla_i^\perp \Psi_j^k (s,s) = \begin{cases}
\pm \sin\omega_{j,1} & i = 1\\
\pm \sin \omega_{j,2} &i = 2,
\end{cases}
\end{align*}
\noindent where $\omega_{j,1}$ and $\omega_{j,2}$ are the angles made with the initial and final links respectively of $\gamma_j(s)$ and the positively oriented tangent line $T_s \d \Omega$. These quantities are clearly invariant under dilation, which is in fact a conformal transformation and the $\pm$ signs only depend on whether the orbits emanating from $s$ originated from orbits in the $TT, TN, NT$, or $NN$ configuration.

\subsection{Integrating the principal symbol}\label{Integrating the principle symbol}
Recall that $\chi_1$ was a smooth cutoff identically equal to $1$ near $[t_j,T_j]$ as discussed in Section \ref{Oscillatory integral representation}. Adding up the various contributions from orbits in $\mathcal{G}_j$, we find that locally near $[t_j,T_j]$, modulo Maslov factors and lower order terms, we have
\begin{align}\label{D1}
\begin{split}
&\int_{0^+}^1 \chi_1(t/r) \d_r \text{Tr} \cos t\sqrt{-\Delta_r} dr =\\
\sum_\pm \frac{\pm 2t(-1)^{j}}{i} \int_{0^+}^1 \int_0^{|\d \Omega_1|} \int_{0}^\infty& \Delta^*(e^{\pm i \tau(t - r \Psi_j)} |\tau|^{3/2} r^{1/2}
\sin \omega_1 \sin \omega_2 |A_j|^{1/2} X\cdot N) d\tau ds_1 dr.
\end{split}
\end{align}
In order to apply the formal manipulation suggested by equation \eqref{intuition hadamard} in Section \ref{Computing the wave trace}, we need to understand how intertwining localization, integration and radial differentiation in \eqref{D1} combine to produce the desired trace formula. In equation \eqref{D1}, we have also assumed that integrating the principal symbol of the radially differentiated wave trace localized near $[rt_j,rT_j]$ in fact gives us the principal symbol of the honest wave trace near $[t_j,T_j]$. This relies on a structural characterization of the subprincipal terms in our parametrix for $S_j$ (see Lemma \ref{indep of t} below). The latter is carried out in Lemma \ref{integrate from delta} of Section \ref{Lower order terms} and we first address the localization problem.

\begin{lemm}\label{integrate from delta}
To obtain the principal symbol of the wave trace localized near lengths in $[t_j,T_j]$ of a given periodic orbit, it suffices to integrate the variational formula localized near $[rt_j,rT_j]$ from $r_0$ to $1$ rather than $0^+$ to $1$, where $r_0 = r_0(j)$ is chosen sufficiently small in terms of $j$. The remainder
\begin{equation}
\int_0^{r_0} \d_r ( \chi_1(t/r) \text{Tr} \cos t \sqrt{-\Delta_r}) dr
\end{equation}
is in fact smooth for $t$ near $[t_j,T_j]$.
\end{lemm}
\begin{proof}
The localized wave trace near $[t_j,T_j]$ is given by
\begin{equation}
\chi_1(t) \text{Tr} \cos t \sqrt{\Delta_1}.
\end{equation}
All geodesic loops of rotation number $1/j$ have length in the interval $[t_j, T_j]$, where $t_j = \inf_{\omega(\gamma) = 1/j} \text{length}(\gamma)$ and $T_j =\sup_{\omega(\gamma) = 1/j} \text{length}(\gamma)$ (recall that $\omega$ is the function which outputs the rotation number of a given periodic orbit). The work in \cite{MaMe82} shows that there exist a $j_0 \in \mathbb{N}$ such that for $j,k > j_0$, the intervals $[t_j, T_j]$ and $[t_k,T_k]$ are disjoint whenever $j \neq k$. Moreover, we also have the asymptotic $T_j - t_j = O(j^{-\infty})$.
As the Euclidean Laplacian scales by $r^{-2}$ under the dilation $x \to rx$, the wave trace also scales by
\begin{align}\label{rescaled trace}
\text{Tr} \cos t\sqrt{-\Delta_{\Omega_r}} = \sum_j \cos t \lambda_j/r,
\end{align}
where $(\lambda_j^2)_j$ are the Dirichlet eigenvalues of $\Delta_{\Omega_1}$. It is clear that the rescaled trace \eqref{rescaled trace} is a smooth family of distributions in the weak topology. Hence, for $0 < r_0 < 1$ to be chosen sufficiently small later, we have
\begin{equation}\label{int localized trace}
\int_{r_0}^{1} \d_r \left( \chi_1(t/r) \text{Tr} \cos t \sqrt{-\Delta_r} \right)dr = \chi_1(t) \text{Tr} \cos t \sqrt{-\Delta_1} - \chi_1(t/r_0) \text{Tr} \cos t \sqrt{-\Delta_{r_0}}.
\end{equation}
Note that $[t_j,T_j] \subset [0,|\d \Omega|)$, so if $t > 2 r_0 |\d \Omega|$, then $\chi_1(t/r_0) = 0$. Hence, we choose $0 < \delta < t_j/(2 |\d \Omega|)$.\\
\\
However, in contrast to \eqref{int localized trace}, the quantity integrated in \eqref{D1} was localized near $[rt_j, rT_j]$ only \textit{after} taking the variation of the wave trace. We see that the integrand in equation \eqref{int localized trace} splits into
\begin{equation}\label{product rule}
\d_r \left(\text{Tr} \cos t \sqrt{-\Delta_r} \chi_j(t/r)\right) = \chi_j(t/r) \d_r \text{Tr} \cos t \sqrt{-\Delta_r} - r^{-2} \chi_j'(t/r) \text{Tr} \cos t \sqrt{-\Delta_r}.
\end{equation}
As $\chi_j$ was identically equal to $1$ on $[t_j,T_j]$, we see that $\chi_j'(\cdot/r)$ vanishes on $[rt_j,rT_j]$. Hence, $\text{Supp} \chi_j'(\cdot/r) \cap \text{SingSupp Tr} \cos \cdot \sqrt{-\Delta_r} = \emptyset$ and integrating the second term on the righthand side of \eqref{product rule} in $r$ contributes only smooth errors in $t$ to the wave trace near $[t_j,T_j]$.
\end{proof}

\begin{rema}
For any Schwartz function $\phi$ on $\mathbb{R}$, we have
\begin{align}\label{action on a test function}
\text{Tr} \cos t \sqrt{-\Delta_{\epsilon}} (\phi) = \sum \int_\mathbb{R} \phi(t) \cos t \lambda_j /r.
\end{align}
Integrating by parts $4k$ times in the formula \eqref{action on a test function}, we have
\begin{align}\label{int by parts}
\text{Tr} \cos t \sqrt{-\Delta_{\epsilon}} (\phi) =  \sum_j \left(\frac{r}{\lambda_j}\right)^{4k} \int_\mathbb{R} \phi^{(4k)}(t) \cos t \lambda_j /r.
\end{align}
For $k$ sufficiently large, \eqref{int by parts} can be seen to be absolutely convergent by Weyl's law and an application of the dominated convergence theorem then shows that
$$
\lim_{\epsilon \to 0^+} \text{Tr} \cos t \sqrt{-\Delta_{\epsilon}} = 0
$$
in the sense of distributions. In the Neumann case, $\lambda_0= 0$, which implies by the same logic that
$$
\lim_{\epsilon \to 0^+} \text{Tr} \cos t \sqrt{-\Delta_{\epsilon}} = 1.
$$
This only contributes a smooth error to the trace formula in \eqref{Main theorem}. However, the proof of Lemma \ref{integrate from delta} above actually resolves the issues near $r = 0$ at the expense of smooth errors. Nonetheless, it seemed conceptually more clear to motivate the variational approach by anticipating that the highly oscillatory terms near $r = 0$ would vanish in the limit.
\end{rema}

\subsection{Proof of Theorem \ref{Main theorem}}
We now make the change of variables $\xi = r \tau$, which turns equation \eqref{D1} into
\begin{align}\label{D2}
\sum_{\pm} \frac{(-1)^j2t}{\pm i} \int_{r_0}^1 \int_0^{|\d \Omega_1|} \int_{0}^\infty \Delta^*(e^{\pm i \xi t/r} e^{-i \xi \Psi_j} |\xi|^{3/2} \sin \omega_{j,1} \sin \omega_{j,2} |A_m^k|^{1/2} X\cdot N) d\xi ds_1 \frac{dr}{r^2}.
\end{align}
Separating out the oscillatory $dr$ integral, we see that
$$
\int_{r_0}^1 e^{\pm i\xi t/ r} r^{-2}dr =  \frac{e^{\pm i\xi t} - e^{\pm i\xi t/r_0}}{\pm i \xi t}.
$$
This simplifies \eqref{D2} to
\begin{align}\label{D3}
\begin{split}
 \sum_{\pm} \pm {2} (-1)^{j+1} \int_0^{|\d \Omega_1|} \int_{0}^\infty \Delta^*( e^{\pm i \xi (t- \Psi_j)} |\xi|^{1/2} \sin \omega_1 \sin \omega_2 |A_j|^{1/2} X \cdot N) d\xi ds\\
\mp  2 (-1)^{j+ 1} \int_0^{|\d \Omega_1|} \int_{0}^\infty \Delta^* (e^{\pm i \xi (t/r_0  - \Psi_j)} |\xi|^{1/2} \sin \omega_1 \sin \omega_2 |A_j|^{1/2} X\cdot N) d\xi ds.
\end{split}
\end{align}
We now apply the nonstationary phase lemma to the second term in \eqref{D3}; the critical set is given by $\{t = r_0 \Psi_j(s,s)\}$ and since all terms are localized by $\chi_1$ near $[t_j,T_j]$, our choice of $r_0$ in Lemma \ref{integrate from delta} shows that $r_0 \Psi_j(s,s) \notin [t_j,T_j]$ for all $s \in \d \Omega$. Hence, this second term is actually smooth. Equation \eqref{D3} is exactly the quantity appearing in Theorem \ref{Main theorem} without Maslov factors. In Section \ref{Maslov factors}, it is explained why the real part of the integral is taken in Theorem \ref{Main theorem}, ultimately owing to conjugacy of the Maslov factors on the two branches of $\Gamma^j$. In Section \ref{Lower order terms} below, it is verified that radially integrating the principal symbol of the variational wave trace indeed gives the principal symbol of the honest wave trace (cf. remarks preceding Lemma \ref{integrate from delta}).

\subsection{Lower order terms}\label{Lower order terms}
\noindent As mentioned in the remarks following equation \eqref{D1}, we have to justify that integrating the principal symbol of $\d_r \text{Tr} \cos t \sqrt{- \Delta_{\Omega_r}}$ in the variable $r$ gives the principal symbol of $\text{Tr} \cos t \sqrt{-\Delta_{\Omega_1}}$. The parametrix
$$
S(t) = \sum_j S_j(t)
$$
described in \cite{Ch76} approximates the wave kernel modulo $C^\infty$ and Theorem \ref{parametrix with sum} gives us an oscillatory integral expansion for each $S_j$ in terms of the operators $S_m^k$ modulo lower order terms. In addition to the principal symbol computed in Lemma \ref{symbol prop}, $S_m^k$ in fact has a complete parametrix on each dilated domain $\Omega_r$:
\begin{align}
S_{m,\pm}^k(t,x,y; r) = \int_{0}^{+\infty} e^{\pm i \tau (t-\Psi_m^k(x,y;r))} a_{m,k,\pm}(t, \tau, x, y;r) d \tau,
\end{align}
where $a_{m,k,\pm} \in S_{\text{cl}}^{1/2}(\mathbb{R} \times \mathbb{R} \Omega \times \Omega)$ is a classical elliptic symbol of order $1/2$, implicitly depending on $r$. The full symbol has an asymptotic expansion in negative half integer powers of $\tau$:
\begin{align}\label{full symbol}
a_{m,k,\pm}(t,\tau, x,y; r) \sim \sum_{p = 0}^\infty a_{m,k,p,\pm}(t,x,y;r) |\tau|^{1/2 - p}.
\end{align}
We already know that $\Psi_j^k(x,y;r) = r \Psi_j^k(x/r,y/r;1)$ is homogeneous of degree one in $r$. The potential problem is that integrating lower order terms $a_{m,k,p,\pm}$ in the radial variable $r$, for $p \geq 1$ could actually produce large powers of $\tau$ which contribute to the principal symbol of $\text{Tr} \cos t \sqrt{-\Delta_1}$. We now show that this is not the case, by taking advantage of certain homogeneities for the dilated kernels and symmetries in the subprincipal terms on the critical set.
\begin{lemm}\label{indep of t}
The full symbol \eqref{full symbol} can be made independent of $t$ for $x, y$ in the open region where the parametrix \ref{parametrix with sum} is valid. Moreover, each term in the asymptotic expansion contributes $O(|\tau|^{-p-1})$ to the full symbol of the wave trace as $\tau \to \pm \infty$.
\end{lemm}
\begin{proof}
On the canonical relation $\Gamma_{\pm}^{m,k}$, we have $t = \Psi_m^k$. Hence, substituting $\widetilde{a}_{m,k,p,\pm}(x,y;r) := a_{m,k,p,\pm}(\Psi_j^k(x,y), x,y;r)$ in \eqref{full symbol} defines the kernel of an operator which differs from $S(t)$ by a smoothing operator microlocally near $\Gamma_{\pm}^{m,k}$. Note also that the wave kernel has the following homogeneity property:
\begin{align}\label{homogeneity}
S(t,x,y;r) = r^{-1} S(t/r, x/r, y/r; 1) \qquad x,y \in \Omega_r.
\end{align}
This can be seen by computing
$$
(\d_t^2 - \Delta_{r}) r S(t/r,x/r,y/r;1) =  r^{-1} (\d_s^2 - \Delta_{1}) S(s,z, w ; r) \Bigg|_{\substack{s = t/r\\ z = x/r\\ w = y/r}} = 0,
$$
for $x,y \in \Omega_r$. As both kernels in \eqref{homogeneity} solve the homogeneous wave equation with the same initial conditions
\begin{align}
\begin{split}
S(0,x,y;r) &= {r}^{-1} S(0,x/r,y/r;1) = 0,\\
\frac{\d S}{\d t} (0,x,y ; r) &= r^{-1} \frac{\d}{\d t} S(t/r,x/r,y/r ; 1 ) \Big|_{t = 0} = \delta (x-y),
\end{split}
\end{align}
uniqueness for the wave equation then implies the asserted homogeneity. Upon equating the phase functions via a change of variables, we obtain the symbolic equation
\begin{equation}\label{homogeneity of symbols}
a_{m,k,p,\pm} (t, rx, ry; r) = r^{1/2- p} a_{m,k,p,\pm}(t/r, x, y;1) = r^{1/2 - p} \widetilde{a}_{m,k,p,\pm}(x,y;1).
\end{equation}
The integrand appearing in the variation of the wave trace \eqref{Hadamard reference} is $\nabla_1^\perp \nabla_2^\perp S_r(t,q,q)$. We should actually replace the  gradients with the correctly scaled ones $r^{-1}\nabla_i^\perp$, which are the differential operators corresponding to unit length vector fields for each $r$. We calculate that
\begin{align}\label{wave gradient}
\begin{split}
&\left(r^{-2} \nabla_1^\perp \nabla_2^\perp S_r \right) (t,rs,rs)\\
& = \sum_{p \geq 0, \pm} \int_0^\infty e^{\pm i\tau (t - r \Psi_m^k)} (\nabla_1^\perp \Psi_m^k)(\nabla_2^\perp \Psi_m^k) a_{m,k,p,\pm}(t,rs,rs;r) |\tau|^{3/2 - p} d\tau.
\end{split}
\end{align}
Changing variables $\xi = r\tau$ and using the homogeneity in equation \eqref{homogeneity} in addition to that of $\nabla \Psi_m^k$, we see that integrating \eqref{wave gradient} in $r$ gives
\begin{align}
\begin{split}
&\int_{r_0}^1 \frac{\d}{\d r} \text{Tr} \cos t\sqrt{-\Delta_r} dr \\
&= \sum_{p \geq 0, \pm} \int_0^{|\d \Omega|} \int_{r_0}^1 \int_0^\infty e^{\pm i \xi (t/r - \Psi_m^k)} (\nabla_1^\perp \Psi_m^k)(\nabla_2^\perp \Psi_m^k) \widetilde{a}_{m,k,p,\pm}(s,s;1) |\xi|^{3/2 - p} \frac{d\xi}{r^2} dr ds.
\end{split}
\end{align}
Proceeding exactly as in equations \eqref{D2} and \eqref{D3} from the previous section, we conclude the proof of the lemma.
\end{proof}
\begin{rema}
While Lemma \ref{indep of t} shows that integrating the subprincipal terms $a_p$ in the radial variable, for $p \geq 1$, does not contribute large powers of $\tau$, it is still possible that smooth remainders blow up near $r = 0$. However, replacing $0^+$ by $r_0$ as in Lemma \ref{integrate from delta} resolves this issue.
\end{rema}
\subsection{Maslov factors}\label{Maslov factors}
To explicitly compute the Maslov factors $\sigma_j^\pm$ on $\Gamma_{\pm}^{j}$, we use an argument due to Keller (\cite{Keller}), following the presentation in \cite{HassellHillairetFord}. The free wave propagator $U(t) = e^{-i t \sqrt{-\Delta}}$ on $\mathbb{R}^2$ has an integral kernel given by
\begin{align}
U(t,x,y) = (2\pi)^{-2} \int_{\mathbb{R}_\xi^2} e^{i(\langle x-y, \xi \rangle - |\xi| t)}d\xi|dx \wedge dy|^{1/2},
\end{align}
considered as a distributional half density (cf. Section \ref{Fourier Integral Operators}). By changing variables and applying the method of stationary phase, it is shown in \cite{HassellHillairetFord} that the principal symbol of $U(t)$ on $N^*\{|x-y| = t\} = \Gamma_{\pm}^0$ is
\begin{align}\label{Mas}
e^{-i\pi/4} \left(\frac{\tau}{t}\right)^{1/2} |d\tau \wedge dx \wedge dy|^{1/2}.
\end{align}
Hence, the Maslov factors on $\Gamma_{\pm}^0$ are given by $\sigma_0^{\pm} = \pm 1$. The arguments in \cite{GuMe79b}, which are in turn based on the construction in \cite{Ch76}, then show that after a reflection at the boundary, the Maslov factors remain unchanged. Hence, $\sigma_j^\pm = \pm 1$ for all $j \in \mathbb{Z}$. Both of the propagators corresponding to $\Gamma_{+}^j$ and $\Gamma_-^j$ contribute to the wave trace singularity near $[t_j,T_j]$, owing to the two branches of propagation:
\begin{align}\label{two branches}
S(t) = \frac{e^{i t \sqrt{-\Delta}} - e^{- i t \sqrt{-\Delta}}}{2 i \sqrt{-\Delta} }.
\end{align}
Hence, we multiply each $\pm$ branch of the Lagrangian distributions $S_{j, \pm}$ by $e^{\pm i \pi/ 4}$, which explains the real parts taken in Theorem \ref{Main theorem}, Theorem \ref{Rational caustic theorem} and Corollary \ref{Ellipse corollary}. Noting that the principal symbol of $\sqrt{-\Delta}$ is $|\xi|$ and $\tau = \pm |\eta| = \pm |\xi|$ on the canonical relation $\Gamma_{\pm}^0$, equations \eqref{Mas} and \eqref{two branches} also show that the principal symbol of $S(t)$ on $\Gamma^0$ has the correct order.

\subsection{Proofs of Theorem \ref{Rational caustic theorem} and Corollary \ref{Ellipse corollary}}
We now show how Theorem \ref{Rational caustic theorem} and Corollary \ref{Ellipse corollary} follow directly from Theorem \ref{Main theorem}. Recall from Section \ref{Billiards} that a caustic is a smooth curve $\mathcal{C}$ in $\Omega$ such that every link tangent to $\mathcal{C}$ remains tangent to $\mathcal{C}$ after a reflection at the boundary. It is well known (see \cite{Katok}, \cite{Tabachnikov}) that if one periodic orbit is tangent to a caustic $\mathcal{C}$, then every orbit tangent to $\mathcal{C}$ is in fact periodic, with the same period (number of bounces) and winding number. In this case, we say $\mathcal{C}$ is a rational caustic and use $\omega(\mathcal{C})$ to denote the (rational) rotation number of any orbit $\gamma$ tangent to $\mathcal{C}$. Hence, rational caustics correspond to highly degenerate periodic orbits in the sense that they are not isolated and $1$ is an eigenvalue of the Poincar\'e map. Let $L_j$ be the length of a periodic orbit which is tangent to a caustic $\mathcal{C}_j$, making $j$ reflections at the boundary and one rotation. As all orbits tangent to $\mathcal{C}$ are periodic orbits, each $q \in \d \Omega$ is a critical point of the loop function, i.e. $\d_q \Psi_j(q,q) = 0$. In this case, the length function $\Psi_j(q,q)$ appearing in the phase of $\text{Tr} \cos t \sqrt{-\Delta}$ (as in Theorem \ref{Main theorem}) is the constant function $L_j$, i.e. every $j$ reflection loop is in fact a periodic orbit of length $L_j$. Assuming the noncoincidence condition \eqref{NCC} on $\Omega$, all periodic orbits of length $L_j$ arise in this way. Tangency to $\mathcal{C}$ also implies that the angles $\omega_1$ and $\omega_2$ in the amplitude are equal. Hence, the wave trace in Theroem \ref{Main theorem} is given  by
\begin{align}\label{almost caustic}
\Re \left\{(-1)^j e^{-i \pi / 4} {4} \int_{\d \Omega} \int_{0}^\infty e^{i \xi (t- L_j)} |\xi|^{1/2} \sin^{3/2} \omega_1 \left|\frac{\d \omega_1}{\d \theta}\right|^{1/2} X \cdot N d\xi dq\right\}.
\end{align}
As the $dq$ and $d\xi$ integrals can be separated, we obtain the Fourier transform of the homogeneous distribution
\begin{align*}
\chi_+^{3/2}(\xi) = \begin{cases}
\xi^{3/2} & \xi \geq 0\\
0 & \xi < 0
\end{cases}
\end{align*}
One can define $\chi_+^a$ similarly as an $L_\text{loc}^1$ function for $\Re a > -1$ and these distributions fact be analytically continued to a larger region of $a \in \mathbb{C}$. It is shown in \cite{Ho90} (Chapter 7) that the Fourier transform of $\chi_+^a$ (with dual variable $t$) is given by $e^{-i\pi (a+1)/2} (t - i 0)^{-a-1}$. The proof of Theorem \ref{Rational caustic theorem} is concluded by evaluating the Fourier transform in equation \eqref{almost caustic} at the point $t - L_j$.
\\
\\
In the case of an ellipse, the billiard flow is known to be completely integrable and each confocal ellipse is in fact a caustic. Moreover, Poncelet's thoerem (see \cite{Poncelet1}, \cite{Poncelet2}) implies that all periodic orbits tangent to a given confocal ellipse have the same length. Hence, all periodic orbits in the ellipse correspond to rational caustics, so Theorem \ref{Rational caustic theorem} applies. Calculations from the author's previous work using action-angle coordinates and Jacobi elliptic function theory allow for the explicit computation of $\d \omega / \d \theta$ appearing in the integrand of \eqref{almost caustic} (see Section 5.5 of \cite{Vig18}). As the boundary of an ellipse
$$
\left\{(x,y): \frac{x^2}{a^2} + \frac{y^2}{b^2} \leq 1 \right\}
$$
is easily parametrized by $(a \cos \phi, b \sin \phi)$ for $\phi \in [0,2 \pi)$, the quantity $X(q) \cdot N(q) dq$ can be explicitly calculated. Combining these observations with Theorem \ref{Rational caustic theorem}, we obtain the formula appearing in Corollary \ref{Ellipse corollary}.

\section{An auxiliary check on the order of $S_j^k$}\label{Check on the order}
\noindent Given a discrepancy in the works \cite{MaMe82} and \cite{Popov1994}, we provide an additional check on the order of $a_0^j$ in Theorem \ref{Main theorem}. Assume $\Omega$ satisfies the noncoincidence condition \eqref{NCC} and let $\rho \in \mathcal{S}(\mathbb{R})$ be a test function such that $\text{Supp} {\rho} \subset [t_j - \epsilon ,T_j + \epsilon]$, where $\epsilon$ is sufficiently small to ensure $\text{Supp} \rho \cap \text{LSP}(\Omega) = [t_j, T_j]$. Let $L_j \in [t_j, T_j]$ denote the length of a periodic orbit of rotation number $1/j$. We will compute the quantity
\begin{align}\label{SP1}
\mathcal{F}({{\rho}(t) \text{Tr} \cos t \sqrt{-\Delta}} ) = \int e^{-i\lambda t} \rho(t) \text{Tr} \cos t \sqrt{-\Delta} dt
\end{align}
in two different ways. Recalling our parametrix for $\cos t \sqrt{-\Delta}$, we see that \eqref{SP1} is given by
\begin{align}\label{SP2}
\sum_{\pm} \int_{\mathbb{R}_t} \int_{0}^\infty \int_{\d \Omega} e^{\pm i \tau (t - \Psi_j(q,q)) - \lambda t} \rho(t) \tau^m a^j(q) dq d \tau dt,
\end{align}
where $m$ is the purported order of $a^j(q)$. Changing variables by $\xi = \tau/\lambda$, we see that \eqref{SP2} becomes
\begin{align}\label{SP5}
\sum_{\pm} \int_{\mathbb{R}_t} \int_{0}^\infty \int_{\d \Omega} e^{\pm i \lambda \xi (t - \Psi_j(q,q)) - \lambda t} \rho(t) \lambda^{m+1} \xi^m a_0(q) dq d \xi dt.
\end{align}
To understand the $\lambda$ asymptotics of this oscillatory integral, we apply the method of stationary phase. First assume $L_j$ is a simple nondegenerate length. On the critical set, $ d_{t,\xi,q} (\xi(t - \Psi_j(q,q)) - t)  = 0$, which implies
\begin{align*}
\begin{cases}
t = \Psi_j(q,q)\\
\xi = 1\\
d_q \psi_j(q,q) = 0.
\end{cases}
\end{align*}
Hence, \eqref{SP1} is given by
\begin{align}\label{SP4}
(2\pi)^{3/2} \lambda^{m - 1/2} \sum_{q: d_q \Psi_j(q,q) = 0} \frac{e^{ i\pi/4\text{sgn Hess}(\Psi_j(q,q)) }\rho(\Psi_j(q,q))a_0(q)}{| \d_q^2(\Psi_j(q,q))|^{1/2}} + o(\lambda^{m-1/2}).
\end{align}
Recall that periodic orbits of rotation number $1/j$ arise from critical points of the $j$-loop function, i.e. $d_q \Psi_j(q,q) = 0$ implies that the geodesic loop of $j$ reflections based at $q$ is actually a periodic orbit. Since $L_j$ was assumed to be simple, corresponding to a unique nondegenerate orbit, there are precisely $j$ such boundary points $q$ and the sum in \eqref{SP4} is actually finite. It is shown in Theorem 3 of \cite{KozTresh89} that nondegeneracy of $\gamma$ also implies $\d_q^2(\Psi_j(q,q)) \neq 0$.
Now, by the formulas in \cite{GuMe79b}, we know that for a simple length $L_j$ corresponding to a nondegenerate periodic orbit of $j$ reflections, the leading asymptotic of the wave trace modulo Maslov factors is given by
\begin{align}\label{GuMePoisson}
\sum \cos t \lambda_j = \frac{L_j}{|\det(I - P_\gamma)|^{1/2}} (t - L_j + i0)^{-1} \mod L_{loc}^1(\mathbb{R}),
\end{align}
where $P_\gamma$ is the Poincare map associated to the unique periodic orbit $\gamma$ of length $L_j$.  Formulas in \cite{Ho90} tell us that the Fourier transform of the righthand side of \eqref{GuMePoisson} is a constant multiple of the Heaviside function. Comparing degrees of homogeneity in \eqref{SP4} and \eqref{GuMePoisson} immediately implies that $m = 1/2$.
\\
\\
If there are infinitely many critical points of $\Psi_j(q,q)$ in the phase of \eqref{SP4}, the analysis is more subtle. For example, Poncelet's theorem for elliptical billiards actually implies that periodic orbits of a fixed length and rotation number come in one parameter families. For $j$ sufficiently large, every boundary point is the base point for a unique periodic orbit tangent to a single confocal ellipse, making $j$ reflections and a single rotation. The lengths of these orbits are independent of the base point, which implies $d_q \Psi_j(q,q)$ vanishes identically. In this case, $\Psi_j(q,q) = L_j$ and applying stationary phase to \eqref{SP5} yields
\begin{align}\label{SP3}
2\pi \lambda^{m} e^{i \lambda L_j}  \rho(L_j) \int_{\d \Omega}  a_0(q)dq  + o\left( \lambda^{m}\right).
\end{align}
While the expansion \eqref{GuMePoisson} is no longer valid for high length spectral multiplicity, the formulas in \cite{HeZe12} and \cite{Vig18} show that the \textit{variation} of the wave trace near such a period in the length spectrum is a distribution of the form $d_j \Re (t- L_j + i0)^{-5/2}$ for some constant $d_j$.
Formally, the wave trace has one higher degree of regularity than its variation, which suggests that the wave trace is of the form $c_j (t - L_j +i 0)^{-3/2}$, in agreement with Theorem \ref{Rational caustic theorem} and Corollary \ref{Ellipse corollary}. Comparing this with the asymptotics in \eqref{SP3} and formulas in \cite{Ho90} for the Fourier transform of homogeneous distributions, we see again that $m = -1 - (-3/2) = 1/2$. The order of $a_0^j$ was also confirmed in the recent paper \cite{HeZe19}

\section{Acknowledgements} \noindent The author would like to thank Hamid Hezari for his many suggestions throughout this project.

\nocite{Arnold}
\nocite{Zw12}
\nocite{Ho485}
\nocite{PeSt17}

\bibliographystyle{alpha}

redpajama_set_name: RedPajamaArXiv