Robin Spectral Rigidity of the Ellipse

TL;DR

This paper proves that smooth, symmetry-preserving Robin isospectral deformations of an ellipse are trivial, showing such deformations must be flat at the initial state and do not exist analytically.

Contribution

It introduces a novel approach combining Hadamard's formula and wave trace analysis to establish rigidity of Robin spectral data for ellipses.

Findings

No non-trivial smooth isospectral deformations preserve symmetry.

Deformations must be flat at the initial ellipse configuration.

Constructs explicit wave propagator parametrix near rotation orbits.

Abstract

In this paper, we investigate $C^1$ isospectral deformations of the ellipse with Robin boundary conditions, allowing both the Robin function and domain to deform simultaneously. We prove that if the deformations preserve the reflectional symmetries of the ellipse, then the first variation of both the domain and Robin function must vanish. Reparameterizing allows us to show that such smooth deformations must be flat at $\epsilon = 0$. In particular, there exist no such analytic isospectral deformations. The key ingredients are a version of Hadamard's variational formula for variable Robin boundary conditions and an oscillatory integral representation of the wave trace variation which uses action angle coordinates for the billiard map. For the latter, we in fact construct an explicit parametrix for the wave propagator in the interior, microlocally near orbits of rotation number $1/j$.

Equations

$$\{\begin{array}{ll}-\Delta u_{\epsilon }={\lambda }^2(\epsilon )u_{\epsilon },& x\in {\Omega }_{\epsilon },\\ \frac{\partial u_{\epsilon }}{\partial \nu }=K_{\epsilon }{u}_{\epsilon },& x\in \partial {\Omega }_{\epsilon }.\end{array}$$

$${\int }_{\partial {\Omega }_0}-⟨{\nabla }_2^T{G}_0(\lambda ,x,q),{\nabla }_1^T{G}_0(\lambda ,q,y)⟩\dot{\rho }+({\lambda }^2\dot{\rho }$$

$$\delta \text{Tr}\cos t\sqrt{-{\Delta }_{\epsilon }}$$

$$c_j={\int }_0^{2\pi }P({\lambda }_j)Q(\phi )\frac{\dot{\rho }(\phi )d\phi }{\sqrt{(b^2+(a^2-b^2){\sin }^2\phi )-{\lambda }_j^2}}.$$

$$\frac{x^2}{a^2-{\lambda }_j^2}+\frac{y^2}{b^2-{\lambda }_j^2}=1,$$

$$\delta \text{Tr}\cos t\sqrt{-{\Delta }_{\epsilon }}=\widetilde{c_j}t\mathrm{Re}\left\{e^{\sigma i\pi /4}(t-T-i{0}^{+}{)}^{-1/2}\right\}+L.O.T.,$$

$$\widetilde{c_j}={\int }_0^{2\pi }\tilde{P}({\lambda }_j)\tilde{Q}(\phi )\frac{\dot{K}(\phi )d\phi }{\sqrt{(b^2+(a^2-b^2){\sin }^2\phi )-{\lambda }_j^2}}.$$

$${\int }_{{\Omega }_{\epsilon }}⟨\nabla u_{\epsilon },\nabla {\phi }_{\epsilon }^{-1\ast }v⟩-{\lambda }_j^2{u}_{\epsilon }{\phi }_{\epsilon }^{-1\ast }vdV={\int }_{{\Omega }_{\epsilon }}{f}_{\epsilon }{\phi }_{\epsilon }^{-1\ast }vdV+{\int }_{\partial {\Omega }_{\epsilon }}{K}_{\epsilon }{u}_{\epsilon }{\phi }_{\epsilon }^{-1\ast }vd{q}_{\epsilon },$$

$$e(u,v)=⟨\nabla u,\nabla v⟩-{\lambda }^{2}uv$$

$$\delta u_{\varepsilon}=\frac{d}{d\varepsilon}\Big{|}_{\varepsilon=0}u_{\varepsilon}.$$

$$\delta u(\alpha)=\frac{d}{d\varepsilon}\Big{|}_{\varepsilon=0}u_{\varepsilon}(\alpha),$$

$$\theta_{X}u_{\varepsilon}=\frac{d}{d\varepsilon}\Big{|}_{\varepsilon=0}\varphi_{\varepsilon}^{*}u_{\varepsilon}.$$

$${\theta }_{X}u=\delta u_{\epsilon }+X{u}_0.$$

$${\phi }_{\epsilon }^{\ast }{u}_{\epsilon }-u_0={\phi }_{\epsilon }^{\ast }(u_{\epsilon }-u_0)+({\phi }_{\epsilon }^{\ast }-1)u_0,$$

$${\int }_{{\Omega }_0}{\phi }_{\epsilon }^{\ast }e(u_{\epsilon },{{\phi }_{\epsilon }^{-1}}^{\ast }v){\phi }_{\epsilon }^{\ast }dV={\int }_{{\Omega }_0}{\phi }_{\epsilon }^{\ast }{f}_{\epsilon }v{\phi }_{\epsilon }^{\ast }dV+{\int }_{\partial {\Omega }_{\epsilon }}{K}_{\epsilon }{u}_{\epsilon }{{\phi }_{\epsilon }^{-1}}^{\ast }vd{q}_{\epsilon }.$$

$${\int }_{{\Omega }_0}{e}_{\epsilon }({\phi }_{\epsilon }^{\ast }{u}_{\epsilon },v){\phi }_{\epsilon }^{\ast }dV={\int }_{{\Omega }_0}{\phi }_{\epsilon }^{\ast }{f}_{\epsilon }vdV+{\int }_{\partial {\Omega }_{\epsilon }}{K}_{\epsilon }{u}_{\epsilon }{{\phi }_{\epsilon }^{-1}}^{\ast }vd{q}_{\epsilon },$$

$$e_{\epsilon }(u,v)={\phi }_{\epsilon }^{\ast }(e({{\phi }_{\epsilon }^{-1}}^{\ast }u,{{\phi }_{\epsilon }^{-1}}^{\ast }v))$$

$$e_{\epsilon }=\sum _{|\alpha |,|\beta |\le 2}{c}_{\alpha ,\beta }(\epsilon ,x,y)D_x^{\alpha }{D}_y^{\beta },$$

$$\begin{array}{rl}{\int }_{{\Omega }_0}& (e({\theta }_{X}u,v)+({\theta }_{X}e)(u_0,v))dV+e(u_0,v){\theta }_{X}dV\\ & ={\int }_{{\Omega }_0}v{\theta }_{X}fdV+vf{\theta }_{X}dV+\delta {\int }_{\partial {\Omega }_{\epsilon }}{K}_{\epsilon }{u}_{\epsilon }{\phi }_{\epsilon }^{-1\ast }vd{q}_{\epsilon },\end{array}$$

$$\theta_{X}e=\frac{d}{d\varepsilon}\Big{|}_{\varepsilon=0}e_{\varepsilon}$$

$$({\theta }_{X}e)(u,v)=Xe(u,v)-e(Xu,v)-e(u,Xv).$$

$$e_{\epsilon }(u,v)={\phi }_{\epsilon }^{\ast }{w}_{\epsilon },$$

$$\frac{d}{d\varepsilon}\Big{|}_{\varepsilon=0}\varphi_{\varepsilon}^{*}w_{\varepsilon}=Xw_{0}+\delta w_{\varepsilon}.$$

$$\delta w_{\epsilon }$$

$$=⟨\nabla (-Xu),\nabla v⟩+⟨\nabla u,\nabla (-Xv)⟩-{\lambda }^2(-Xu)v-{\lambda }^{2}u(-Xv)$$

$$=-e(Xu,v)-e(u,Xv).$$

$$\begin{array}{rl}\delta {\int }_{{\Omega }_0}{e}_{\epsilon }({\phi }_{\epsilon }^{\ast }{u}_{\epsilon },v){\phi }_{\epsilon }^{\ast }dV& ={\int }_{{\Omega }_0}e(\dot{u},v)-e(u_0,Xv)dV+{\int }_{\partial {\Omega }_0}e(u_0,v)X_{\nu }d{q}_0\\ & ={\int }_{{\Omega }_0}v\dot{f}-(Xv)fdV+{\int }_{\partial {\Omega }_0}vf{X}_{\nu }d{q}_0+\delta {\int }_{\partial {\Omega }_{\epsilon }}{K}_{\epsilon }{u}_{\epsilon }{{\phi }_{\epsilon }^{-1}}^{\ast }vd{q}_{\epsilon }.\end{array}$$

$${\int }_{\partial {\Omega }_{\epsilon }}{K}_{\epsilon }{u}_{\epsilon }{\phi }_{\epsilon }^{-1\ast }vd{q}_{\epsilon }={\int }_{\partial {\Omega }_0}{K}_{\epsilon }(x+{\rho }_{\epsilon }(x){\nu }_x)u_{\epsilon }(x+{\rho }_{\epsilon }(x){\nu }_x)v(x)d{q}_{\epsilon }.$$

$$\delta d{q}_{\epsilon }=\kappa \dot{\rho }d{q}_0,$$

$$\frac{d}{d\epsilon }$$

$${\int }_{{\Omega }_0}$$

Full text

1. Introduction

In this paper, we prove infinitesimal spectral rigidity of the ellipse with $C^{1}$ deformations in both the domain and Robin boundary conditions which preserve the symmetries of the ellipse. This means that the first variations of both the domain and Robin function vanish. To make this precise, we consider the eigenvalue problem for the Laplacian. Let $\Omega_{0}$ be an ellipse and $\varphi_{\varepsilon}$ a $C^{1}$ family of diffeomorphisms defined in a neighborhood of $\Omega_{0}$ for $0\leq\varepsilon<\varepsilon_{0}$, such that $\varphi_{0}=\text{Id}$. Denote $\Omega_{\varepsilon}=\varphi_{\varepsilon}(\Omega_{0})$ and let $K_{\varepsilon}$ be a $C^{1}$ family of smooth functions on $\partial\Omega_{\varepsilon}=\varphi_{\varepsilon}(\partial\Omega_{0})$. The PDE we are interested in is

[TABLE]

The equation (1) is said to have Robin boundary condition with Robin function $K_{\varepsilon}$. As $K_{\varepsilon}\in C^{\infty}(\partial\Omega_{\varepsilon})$, it is shown in [Tay11] that $-\Delta$ is self adjoint on $L^{2}(\Omega_{\varepsilon})$ with densely defined domain $D=\{u\in H^{2}(\Omega_{\varepsilon}):\partial_{\nu}u_{\varepsilon}=K_{\varepsilon}u_{\varepsilon}\}$. The analytic Fredholm theorem guarantees that (1) has nontrivial solutions only for a discrete collection of eigenvalues $\lambda_{j}^{2}\in\mathbb{R}$, comprising the entire spectrum which we denote by $\text{Spec}(\Omega_{\varepsilon})$. We may assume $\varphi_{\varepsilon}(x)=x+\rho_{\varepsilon}(x)\nu_{x}$ for $x\in\partial\Omega_{0}$, where $\nu_{x}$ is the outward unit normal to $\partial\Omega_{0}$ and $\rho_{\varepsilon}$ is a smooth function on $\partial\Omega_{0}$. We also impose the restriction that $\rho_{\varepsilon}$ and $K_{\varepsilon}$ should be invariant under the symmetry group of the ellipse, which is generated by reflections through the coordinate axes and is isomorphic to the Klein four-group. We denote the Laplace operator with boundary conditions above by $\Delta_{\varepsilon}$ and prove the following result:

Theorem 1.1**.**

Let $\varphi_{\varepsilon}(\Omega_{0})=\Omega_{\varepsilon}$ be a $C^{1}$ deformation of the ellipse through smooth domains with $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetries and $K_{\varepsilon}$ a $C^{1}$ family of Robin functions with the same symmetries. If the deformation is isospectral, i.e. $\text{Spec}(\Delta_{\varepsilon})=\text{Spec}(\Delta_{0})$ for $0\leq\varepsilon<\varepsilon_{0}$, then $\dot{\rho}=\dot{K}=0$.

Here, we have written \dot{\rho}(x)=\frac{d}{d\varepsilon}\big{|}_{\varepsilon=0}\rho_{\varepsilon}(x) and \dot{K}(x)=\frac{d}{d\varepsilon}\big{|}_{\varepsilon=0}K_{\varepsilon}(x+\rho_{\varepsilon}(x)\nu_{x}) to denote the first variations, which are studied in more detail in Section 3. We also use the notation \delta=\frac{d}{d\varepsilon}\big{|}_{\varepsilon=0}. We then say the ellipse is infinitesimally spectrally rigid through domains and Robin boundary conditions with the symmetries of an ellipse. The prefix “infinitesimal” means that only the first variation vanishes. If it could be shown that no nontrivial such isospectral deformations exist, we would say the ellipse is spectrally rigid amongst such domains. As a quick corollary to Theorem 1.1, we have:

Corollary 1.2**.**

There are no nontrivial analytic isospectral deformations of the ellipse through $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetric domains and Robin functions.

This follows from a simple reparametrization argument which can be found in Section 3.2 of [HZ12]. One of the main ingredients in our proof is a version of Hadamard’s variational formula for variable Robin boundary conditions, which also appears to be new in the literature:

Theorem 1.3**.**

Let $G_{\varepsilon}$ be the Green’s kernel for the eigenvalue problem with Robin boundary conditions on $\Omega_{\varepsilon}$. Then, for $x,y\in\text{int}(\Omega_{0})$, $\delta G_{\varepsilon}(\lambda,x,y)$ is the distribution

[TABLE]

Here, $\nabla_{i}^{T}$ denotes the tangential derivative in the $i$th spatial variable and $dq$ is the natural line element on $\partial\Omega_{0}$ inherited from the flat metric on $\mathbb{R}^{2}$. The Green’s kernel or Green’s function is the Schwartz kernel of the resolvent $(-\Delta-\lambda^{2})^{-1}$ for $\operatorname{Im}\lambda^{2}>0$. In Section 3, this result is extended to variational formulas for both the even wave trace $\text{Tr}\cos(t\sqrt{-\Delta_{\varepsilon}})$ and simple eigenvalues. In particular, we prove the following in Section 6:

Theorem 1.4**.**

For the ellipse $\Omega_{0}=\{(x,y):\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}\leq 1\}$, with boundary parametrized by $(x,y)=(a\cos\varphi,b\sin\varphi)$, the variation of the wave trace near a simple length $T_{j}\in\text{Lsp}(\Omega_{0})$ corresponding to a caustic of rotation number $1/j$ is given by

[TABLE]

and the constants $c_{j}$ are given by

[TABLE]

Here, L.O.T. denotes lower order distributional terms and $\lambda_{j}$ is the parameter of the confocal ellipse

[TABLE]

to which periodic orbits of length $T_{j}$ are tangent. $P$ and $Q$ are certain nonzero analytic functions. Moreover, if $\dot{\rho}=0$, then the leading order term in the singularity expansion near $T_{j}$ becomes

[TABLE]

where the constants $\widetilde{c_{j}}$ are given by

[TABLE]

$\widetilde{P}$* and $\widetilde{Q}$ are again certain nonzero analytic functions.*

Remark 1.5*.*

In [GM79a], it is shown that there exists a $j_{0}$ such that for all $j\geq j_{0}$, the lengths $T_{j}$ of periodic orbits with rotation numbers $1/j$ are simple. In this context, simplicity of the length $T_{j}$ means that all periodic orbits of length $T_{j}$ are tangent to a single confocal conic section. These orbits have rotation number $1/j$, which means that they make $j$ reflections at the boundary and have winding number $1$. In Section 6, we will send $j\to\infty$ and use the constants $c_{j}$ and $\widetilde{c_{j}}$ to show that $\dot{\rho}$ and $\dot{K}$ vanish.

1.1. Schematic outline

In Section 2, we review the relevant background on the inverse spectral problem for (1). Section 3 is then devoted to proving a version of Hadamard’s variational formula with Robin boundary conditions (Theorem 1.3) and extending this to variational formulas for both the localized wave trace near lengths of periodic orbits (Theorem 3.7) and simple eigenvalues (Theorem 3.8). The billiard map and its special properties in the case of an ellipse are introduced in Section 4. Relevant background information on dynamical systems and various notations are also introduced there which will be used throughout the paper. Section 5 contains the most difficult material. It begins by reviewing Fourier integral operators and Chazarain’s parametrix and then develops an explicit oscillatory integral parametrix for the wave propagator (Theorem 5.7) which is used to prove a variational trace formula (Theorem 5.11) in Section 5.3. An important dynamical lemma (Lemma 5.3) on the structure of approximate geodesic loops is also stated in this section, but its proof is relegated to Section 7. Section 5.5 analyzes the singularity coefficients of the variational trace formula in Theorem 5.11 by using Jacobi elliptic function theory. Section 6 converts the variational trace formula into an elliptic integral and then uses a method of Guillemin and Melrose to show that the first variation of both the Robin function and domain vanish, which completes the proof of Theorem 1.1. Section 7 contains the proof of Lemma 5.3 and is broken up into several intermediate lemmas.

2. Background

The inverse spectral problem has a long history, dating back to Kac in 1967, who asked the famous question “can one hear the shape of a drum?” In [Kac66], a positive answer to this question was obtained in the special case of a Euclidean ball by using heat invariants and the isoperimetric inequality. Almost immediately afterwards, John Milnor found an example of 16 dimensional isospectral tori which were not isometric. Sunada later generalized this example by using an algebraic method, but the question for planar domains remained open until 1992, when distinct polygonal domains in $\mathbb{R}^{2}$ were found to be isospectral in [GWW92]. The question is still widely open for convex and/or smooth domains, although there has been significant progress. For example, it is proved in [Zel09] that analytic domains with a single isometric involution are spectrally determined assuming some additional generic dynamical constraints on the length spectrum. Melrose also showed that the set of isospectral planar domains is precompact in the $C^{\infty}$ topology by a careful analysis of the heat invariants (see [Mel07]). This was improved in [OPS88a], [OPS88b] and [OPS88c], where the authors proved genuine $C^{\infty}$ compactness of the isospectral set. This result is based on the Polyakov formula for $\zeta$-regularized determinants and applies to both bounded planar domains and closed surfaces. Recently, Hezari and Zelditch showed in [HZ19] that ellipses of small eccentricity are spectrally determined amongst all smooth planar domains, which is the first positive result in such generality since Kac’s original paper [Kac66]. Thorough surveys of the inverse spectral problem are contained in [Zel14], [Zel04], [DH13] and [Mel96].

Dual to the Laplace spectrum is the so called length spectrum, which is a discrete set of numbers containing the lengths of periodic orbits for the geodesic or billiard flow. The same inverse problem exists: can one determine a manifold up to isometry from its length spectrum? The answer is unfortunately negative, as was shown for the case of constant negative curvature in [Vig80]. However, it is conjectured by Katok and Burns that the marked length spectrum does determine a smooth closed manifold up to isometry ([BK85]). Here, the marked length spectrum also encodes the homotopy classes of periodic geodesics. Marked length spectral rigidity was recently shown in [GL18] for Anosov manifolds. The relationship with the Laplace spectrum is contained in the Poisson relation, which tells us that the singularities of the wave trace are a subset of the length spectrum. Assymptotic formulas near the singularities are given by the Selberg trace formula for hyperbolic surfaces ([Sel56]), the Duistermaat-Guillemin trace theorem ([DG75]) for general manifolds under a dynamical nondegeneracy condition, and a Poisson summation formula for strictly convex bounded planar domains due to Guillemin and Melrose ([GM79b]). However, since these trace formulae 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. Hence, without length spectral simplicity, there is no way to deduce Laplace spectral information from the length spectrum alone.

While some results on spectral rigidity are known in the chaotic regime (see for example [GK80a], [GK80b] and [PSU14]), very little is known about the completely integrable setting, in which the flow has a maximal number of conserved quantities (see Section 4). In the theory of dynamical systems, a famous conjecture of Birkhoff is that the only strictly convex planar domains with completely integrable billiards are ellipses. While this remains an open conjecture, much progress has been made in the local setting. It is shown in [KS18] that if a rationally integrable billiard table is sufficiently close to an ellipse, then it must be an ellipse. Rational integrability means that for each integer $q\geq 1$, the billiard map has invariant curves of rotation number $1/q$, consisting entirely of periodic points. Using Aubry-Mather theory, the authors then show that ellipses are length spectrally rigid (Corollary 14 of [KS18]).

The ellipse is smooth, convex and has completely integrable dynamics, which makes it an interesting object to study in the context of spectral theory. In fact:

Conjecture 2.1** (Melrose, [Mel96]).**

Ellipses are spectrally determined.

In [dSKW17], the authors show that convex domains with $\mathbb{Z}_{2}$ axial symmetry which are sufficiently close to a circle in $C^{k}$ are length spectrally rigid. These domains include ellipses of small eccentricity. It is shown in [PS92] that generically, convex domains have simple length spectrum and nondegenerate Poincaré map. Combining the results in [dSKW17] and [PS92], it then follows that $\mathbb{Z}_{2}$ symmetric domains close to a circle are $\Delta$ spectrally rigid amongst a generic class of symmetric domains. In [Hez17], these results are extended to the Robin Laplacian, where the Robin function on the boundary is also allowed to deform through smooth functions with the same same $\mathbb{Z}_{2}$ symmetry. Our problem is similar in nature but considers ellipses of arbitrary eccentricity, which might not be close to a circle.

The present article is inspired by [HZ12], [GM79a], [GM79b] and [Pee80]. Guillemin and Melrose proved a version of the Poisson summation formula for bounded planar domains and then used this result in a subsequent article to show that for a fixed ellipse, a $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetric Robin function on the boundary is completely determined by the spectrum of the associated Laplacian. Hezari and Zelditch then proved infinitesimal spectral rigidity for Dirichlet/Neumann boundary conditions, while only letting the domain deform. In their proof, the authors used the symbol calculus in [DG75] to compute the trace of the wave kernel near periodic transversal reflecting rays. Our problem allows both the domain and Robin function to deform simultaneously, which doesn’t allow us to directly employ the results in [GM79a] or [HZ12].

The idea of the proof of Theorem 1.1 is that if the deformation is isospectral, then the variation of the wave trace should also be zero. The Poisson relation, in this case due to Guillemin and Melrose ([GM79a]), tells us that the singularities of the wave trace are contained in the length spectrum $\text{Lsp}(\Omega_{0})$, the set of lengths of periodic trajectories for the broken bicharacteristic (billiard) flow. Using microlocal analysis and in particular, Chazarain’s parametrix, we can localize the wave kernel near the periodic transversal reflecting rays. From this, we obtain a singularity expansion for the wave trace variation, a Fourier integral distribution, near the length spectrum. We then cook up an oscillatory integral which microlocally approximates this distribution by using a special phase function associated to the billiard map. To do this, we actually construct an explicit parametrix for the microlocalized wave kernel near all orbits tangent to a confocal ellipse of rotation number $1/j$ (Theorem 5.7). In particular, this involves finding all orbits making approximately one rotation with a prescribed number of reflections which connect two points in an interior neighborhood of the diagonal of the boundary (Lemma 5.3). This is of independent interest in the theory of dynamical billiards.

In the special case of the ellipse, we can incoorporate action angle coordinates for the billiard map, which allows us to convert the singularity expansion for the wave trace near $\text{Lsp}(\Omega_{0})$ into the product of a nonzero distribution and an elliptic integral as in Theorem 1.4. Since this expansion is valid near any simple length, we can take a special sequence of caustics creeping closer and closer to the boundary. Following the ideas in [GM79a], we send $j\to\infty$ and analyze the coefficients $c_{j}$ and $\widetilde{c_{j}}$ in Theorem 1.4. These coefficients are analytic in the paramater $\lambda^{2}$ and since $\lambda_{j}^{2}\to 0$ as $j\to\infty$, we see that they are actually flat at $\lambda=0$. An application of the Stone-Weierstrass theorem then shows that $\dot{\rho}=0$. Upon substituting $\dot{\rho}=0$, we obtain a new singularity expansion for the subprincipal term with $\dot{K}$ only. The same tricks show $\dot{K}=0$.

3. Variation of the Wave Trace

In this section, we derive variational formulas for the Green’s function, simple eigenvalues and wave trace. The PDE (1) has the weak formulation

[TABLE]

for any $u_{\varepsilon}\in H^{2}(\Omega_{\varepsilon})$ such that $\partial_{\nu}u_{\varepsilon}=K_{\varepsilon}u_{\varepsilon}$ (i.e. in the domain of $-\Delta$ with Robin boundary conditions), $v\in C^{\infty}(\Omega_{0})$ and $f_{\varepsilon}\in\mathcal{E}^{\prime}(\text{int}(\Omega_{\varepsilon}))$ which is an inhomogeneous term for the PDE (1). Here, $dV=dx^{1}\wedge dx^{2}$ is the volume form on $\mathbb{R}^{2}$ and $dq_{\varepsilon}$ is the natural surface measure on $\partial\Omega_{\varepsilon}$ induced from the Euclidian metric. We refer to the the quantity

[TABLE]

as the energy density.

3.1. Variational derivatives

Some care is needed to differentiate the expressions above. We begin by making precise our notion of first variation, following closely the presentation in [Pee80].

Definition 3.1**.**

If $u_{\varepsilon}\in C^{1}([0,\varepsilon_{0}],\mathcal{D}^{\prime}(\Omega_{\varepsilon}))$ is a $C^{1}$ family of distributions, we write $\delta u$, $\delta u_{\varepsilon}$ or $\dot{u}$ for the first variation of $u_{\varepsilon}$ at $\varepsilon=0$, as a distribution in $\Omega_{0}$:

[TABLE]

To simplify notation, for a single function $u$, we will oftentimes write $\delta u=\dot{u}$ and reserve the use of $\delta$ for preceding long formulas. If $\alpha\in C_{0}^{\infty}(\text{int}{\Omega_{0}})$ is a test function, then $\alpha\in C_{0}^{\infty}(\text{int}{\Omega_{\varepsilon}})$ for $\varepsilon\ll 1$ and we can define $\delta u$ by

[TABLE]

i.e. the derivative of a function from $[0,\varepsilon_{0})$ to $\mathbb{R}$. The issue with this definition is that if the supports of the distributions $u_{\varepsilon}$ actually intersect $\partial\Omega_{\varepsilon}$, then the formula above only defines $\delta u_{\varepsilon}$ in the interior of $\Omega_{0}$ and not on the boundary even when $u_{\varepsilon}$ is defined there. For instance, the Green’s kernel is supported near the boundary in the setting of Robin boundary conditions.

To resolve this issue, we follow the ideas in [Pee80], where some geometric heuristics motivate several precise definitions. The set of smooth domains in $\mathbb{R}^{2}$ is an infinite dimensional manifold $\Lambda$, on which the Lie group $\text{Diff}(\mathbb{R}^{2})$ acts. The Lie algebra of $\text{Diff}(\mathbb{R}^{2})$ is the space of smooth vector fields on $\mathbb{R}^{2}$. An initial domain $\Omega_{0}$ and a curve $\varphi_{\varepsilon}$ in $\text{Diff}(\mathbb{R}^{2})$ generate a curve in $\Lambda$, given by $\Omega_{\varepsilon}=\varphi_{\varepsilon}(\Omega_{0})$. Hence, we can associate elements of the Lie algebra to tangent vectors at $\Omega_{0}$. For any given deformation of $\Omega_{0}$, we have an infitessimal generator X=\frac{d\varphi_{\varepsilon}}{d\varepsilon}\big{|}_{\varepsilon=0}. For a fixed $s\in\mathbb{R}$, we can associate to each $\Omega_{0}\in\Lambda$ the fiber $H^{s}(\Omega_{0})$, which is the $L^{2}$ based Sobolev space of order $s$. This defines a smooth vector bundle over $\Lambda$ on which $\text{Diff}(\mathbb{R}^{2})$ acts via pullback (diffeomorphism invariance of the Sobolev spaces). These heuristics motivate the following definition:

Definition 3.2**.**

For a curve $(u_{\varepsilon},\Omega_{\varepsilon})$ in the above vector bundle, we define the Lie derivative to be

[TABLE]

We sometimes drop the $\varepsilon$ subscript and write $\theta_{X}u_{\varepsilon}=\theta_{X}u$ for simplicity. The advantage of this definition is that it is well defined on the boundary of $\Omega_{0}$ for $s>1$, via the Sobolev embedding theorem $H^{s}\hookrightarrow C^{0}$.

Lemma 3.3**.**

If we suppose that $u_{\varepsilon}\in H^{s+1}(\Omega_{\varepsilon})$ is supported away from the boundary and $\theta_{X}u\in H^{s}(\Omega_{0})$, then $\delta u_{\varepsilon}$ exists in $H^{s}(\Omega_{0})$ and

[TABLE]

Proof.

The lemma follows from writing

[TABLE]

dividing both sides by $\varepsilon$ and sending $\varepsilon\to 0$. ∎

Remark 3.4*.*

The formula in Lemma 3.3 is perfectly valid pointwise whenever $x\in\text{int}{\Omega_{0}}$. Both $\theta_{X}$ and $X$ are well defined operators on distributions supported near the boundary, so by setting $\delta u=\theta_{X}u-Xu$, we obtain an extension of Definition 3.1 for $s>1$. From now on, we use $X$ to denote the differential operator acting on distributions on the fixed domain $\Omega_{0}$ and $\theta_{X}$ to denote the Lie derivative acting on distributions or differential forms on $\Omega_{\varepsilon}$ which may also depend on $\varepsilon$.

3.2. A general variational formula

We now derive a variational formula using the weak formulation (2). To obtain an integral equation on the fixed domain $\Omega_{0}$, we pull back the energy density and apply the change of variables formula to (2):

[TABLE]

The pullback of the surface measure $dq_{\varepsilon}$ in the last term on the right is more complicated, so for the time being, we leave it as is. We can rewrite this equation as

[TABLE]

where

[TABLE]

is the conjugated energy density. While $e_{\varepsilon}$ is a composition of operators, it is still of the form

[TABLE]

for some coefficients $c_{\alpha,\beta}$ depending smoothly on $x,y$ and in a $C^{1}$ manner on $\varepsilon$. This will justify use of the product rule when computing $\varepsilon$ derivatives in Lemma 3.5 below. Differentiating (3) in the parameter $\varepsilon$ and setting $\varepsilon=0$ yields

[TABLE]

where the quantity

[TABLE]

is defined analagously to the formula in Definition 3.2.

Lemma 3.5**.**

For $u,v\in H^{2}(\Omega_{0})$, differentiating the conjugated energy density $e_{\varepsilon}$ yields

[TABLE]

Proof.

Consider the family of distributions $w_{\varepsilon}=e({\varphi_{\varepsilon}^{-1}}^{*}u,{\varphi_{\varepsilon}^{-1}}^{*}v)$. Recalling that

[TABLE]

we see by Lemma 3.3 that

[TABLE]

The first term is precisely $Xe(u,v)$ and the second term $\delta w_{\varepsilon}$ is easily calculated by commuting the $\varepsilon$ and $x$ derivatives:

[TABLE]

Commuting derivatives is always valid in the sense of distributions. In order to also restrict to $\varepsilon=0$, note that $[\frac{d}{d\varepsilon},\nabla_{x}]{\varphi_{\varepsilon}^{-1}}^{*}u=0$ in $L^{2}$ and this quantity only involves second derivatives of $u$, first order derivates of $\varphi_{\varepsilon}^{-1}$ in $\varepsilon$, and $x$ derivatives of $\varphi_{\varepsilon}^{-1}$. Here, $[\cdot,\cdot]$ denotes the commutator of two operators. As $u\in H^{2}$ and $\varphi_{\varepsilon}$ is $C^{1}$ with respect to $\varepsilon$, $[\frac{d}{d\varepsilon},\nabla_{x}]{\varphi_{\varepsilon}^{-1}}^{*}u$ is a continuous family (in $\varepsilon$) of $L_{x}^{2}$ functions. Letting $\varepsilon\to 0$ then implies $[\delta,\nabla]{\varphi_{\varepsilon}^{-1}}^{*}u=0$. ∎

Returning to the variation of (3), the Lie derivative of the volume form in equation (4) gives the divergence of $X$, which we would like to convert to a boundary integral, since $X$ is only defined in a tubular neighborhood of $\partial\Omega_{0}$. Using again the formula $\theta_{X}=\delta+X$ and applying the divergence theorem to equation (4) above gives

[TABLE]

Here, $X_{\nu}=\langle X,\nu\rangle$ is the normal component of $X$. For the last term in (5), we can assume the perturbation is in the normal direction and parametrize the boundary by $\partial{\Omega_{0}}\ni x\mapsto x+\rho_{\varepsilon}(x)\nu_{x}\in\partial\Omega_{\varepsilon}$, so that

[TABLE]

We now recall a basic result from differential geometry:

Lemma 3.6**.**

The variation of surface measure is given by

[TABLE]

where $\kappa$ is the curvature of $\partial\Omega_{0}$.

Proof.

A proof using normal coordinates can be found on page 6 of [CM11]. ∎

Hence, differentiating the entire boundary integral (6), we obtain

[TABLE]

Actually, there are two ways to define $\dot{K}$ since it is currently ambiguous as to how to differentiate in $x$ on the hypersurface $\partial\Omega_{\varepsilon}$. The first way involves extending $K_{\varepsilon}$ radially to a function defined on a tubular neighborhood of $\partial\Omega_{\varepsilon}$, so that we may differentiate in $x$ on an open subset of $\mathbb{R}^{2}$. The second way is to define $\dot{K}_{\varepsilon}(x)=d/d\varepsilon|_{\varepsilon=0}K_{\varepsilon}(x+\rho_{\varepsilon}\nu_{x})$. While these two definitions differ pointwise, the integral formula remains the same and we adopt the second definition as it appears more naturally in the proof.

3.3. Variation of Green’s Kernel

Combining equations (5) and (7), we obtain the variational formula

[TABLE]

Recall that $\varphi_{\varepsilon}(x)=x+\rho_{\varepsilon}(x)\nu_{x}$ for $x\in\partial\Omega_{0}$, so that $X_{\nu}=\dot{\rho}$. Integrating by parts and collecting boundary terms, we see that

[TABLE]

If $u$ and $v$ satisfy the PDE (1) with Robin boundary conditions, we have

[TABLE]

Noting that

[TABLE]

we can get rid of one of the $K_{0}^{2}u_{0}v\dot{\rho}$ terms in (8) in exchange for only using tangential derivatives.

To prove Theorem 1.3, we now fix $x_{0}\in\text{int}(\Omega_{0})$ and denote by $G_{\varepsilon}$ be the Greens function on $\Omega_{\varepsilon}$. Formally, the Green’s function is the Schwartz kernel of the resolvent $(-\Delta_{\varepsilon}-\lambda^{2})^{-1}$. Setting $u_{\varepsilon}(x)=G_{\varepsilon}(\lambda,x,y)$, $v(x)=G_{0}(\lambda,x_{0},x)$ and $f_{\varepsilon}(x)=\delta_{x_{0}}(x)$, we obtain

[TABLE]

which is precisely Theorem 1.3 with $x$ replaced by $x_{0}$. As $(\Delta_{\varepsilon}-\lambda^{2})\in\Psi^{2}(\Omega_{\varepsilon})$ is elliptic for $\operatorname{Im}\lambda$ positive, $G_{\varepsilon}$ is a Lagrangian distribution with principal symbol in $S_{1,0}^{-2}$. Hence, $G_{\varepsilon}$ is a family of distributions and $\delta G_{\varepsilon}$ in particular, is a distribution of order $-1$ with wavefront set conormal to the diagonal $\{x=y\}$. As $x,y\in\text{int}(\Omega_{0})$ and $q\in\partial\Omega_{0}$ in Theorem 1.3, the points $(x,q)$ and $(q,y)$ are away from the diagonal, where the distribution is smooth. Hence, the tangential derivatives do not affect the smoothness or integrability. The distribution $\delta G_{\varepsilon}$ can actually be extended up to the boundary using the method of layer potentials, although this is not needed in the remainder of the paper.

3.4. The wave trace and eigenvalues

We now want to find a formula for the variation of the distributional trace of the even wave propagator, $\cos(t\sqrt{-\Delta_{\varepsilon}})$, in terms of that of the Green’s kernel. Recall that the wave propagator $e^{it\sqrt{-\Delta_{\varepsilon}}}$ has a distributional trace in the sense that

[TABLE]

is trace class for any Schwartz function $\varphi$. Its trace is

[TABLE]

where $(\lambda_{j}^{2}(\varepsilon))_{j=1}^{\infty}$ are the eigenvalues of $\Delta_{\varepsilon}$. The sum in (9) can be seen to be convergent via integration by parts combined with Weyl’s law on the asymptotic distribution of eigenvalues. While Weyl’s law is usually stated for Dirichlet or Neumann eigenvalues, the Robin and Neumman asymptotics actually agree up to leading order due to the fact that the boundary operators $\partial_{\nu}-K$ and $\partial_{\nu}$ of the Robin and Neumann Laplacians respectively have the same principal symbol. For example, see [Zay04] or [Ivr16]. Taking real and imaginary parts, the analagous trace formulas hold for the even and odd wave kernels, which we denote by

[TABLE]

respectively. The subscript $R$ here is to denote the Robin boundary conditions. In this section, we shall prove:

Theorem 3.7**.**

The variation of the even wave trace is

[TABLE]

where we have defined

[TABLE]

Here, $\nabla_{i}^{T}$ is again the tangential derivative in the $i$th spacial variable and $\Delta_{2}$ is the Euclidean Laplacian in the second spacial variable. The kernels are first differentiated in the interior using an extension of the tangential vector field and then restricted to the diagonal of the boundary. We will also prove:

Theorem 3.8**.**

If $\lambda_{j}^{2}(0)$ is a simple eigenvalue associated to the $L^{2}$ normalized eigenfunction $\Psi_{j}$, then

[TABLE]

Both theorems are proved together by the same method:

Proof.

Our derivation of the wave trace variation is based on Kato’s variational formulas for sums of eigenvalues in [Kat95]. One has to be careful, as an eigenvalue of higher multiplicity may not be $C^{1}$ in $\varepsilon$. Such eigenvalues can break off to become many different eigenvalues under deformation. However, if we denote by $m(\lambda_{j}^{2})$ the multiplicity of $\lambda_{j}^{2}=\lambda_{j}^{2}(0)$, we will see that the sum $\sum_{1}^{m(\lambda_{j}^{2})}\lambda_{j,k}^{2}(\varepsilon)$ is in fact $C^{1}$ in $\varepsilon$. We actually prove a more general theorem: let $g$ be holomorphic in a neighborhood of the eigenvalue $\lambda^{2}_{j}$ and denote the resolvent operator by $\widetilde{R}_{\varepsilon}(z)=(-\Delta_{\varepsilon}-z)^{-1}$ for $z\notin\text{Spec}(-\Delta_{\varepsilon})$, with Schwartz kernel $\widetilde{G}_{\varepsilon}(z,x,y)$. We write $\widetilde{R}_{\varepsilon}$ and $\widetilde{G}_{\varepsilon}$ since the spectral parameter is $z$ instead of $z^{2}$. Then, by the Cauchy integral formula, we have

[TABLE]

where

[TABLE]

for $\gamma$ a small, positively oriented circle enclosing only the $\lambda_{j}^{2}$ eigenvalue. When $\varepsilon=0$, $T_{g,0}$ is $g$ of the orthogonal projection onto the eigenspace of $\lambda_{j}^{2}$. As the eigenvalues do vary continuously in $\varepsilon$, for $\varepsilon\ll 1$, $T_{g,\varepsilon}$ is the total projector, i.e. $g$ composed with the projection onto the direct sum of the eigenspaces of $\lambda_{j,k}^{2}(\varepsilon)$ for $1\leq k\leq m(\lambda_{j}^{2})$. $T_{g,\varepsilon}$ is in fact a $C^{1}$ family of operators in $\varepsilon$ since the resolvent is. The trace can also be obtained by integrating $\widetilde{G}_{\varepsilon}$ over the diagonal, which combined with equation (10), gives

[TABLE]

As in Section 3.2, we can pull back $\widetilde{G}_{\varepsilon}(z,x,x)dx$ via $\varphi_{\varepsilon}^{*}$ to obtain an integral over a fixed domain, which we then differentiate. We have

[TABLE]

where the last line follows from the definition of $\theta_{X}$ and the standard fact from Riemannian geometry that the Lie derivative with respect to $X$ of the volume form gives the divergence of $X$ times the volume form, i.e. $\theta_{X}(dx)=\text{div}Xdx$. Using the formula for $\theta_{X}$ in Lemma 3.3 and the divergence theorem, we obtain

[TABLE]

We now plug in our variational formula for the Green’s kernel from Theorem 1.3 to see that

[TABLE]

Denote these three integrals by $I_{1},I_{2}$ and $I_{3}$ and let $(\Psi_{j,k})_{k=1}^{m(\lambda_{j}^{2})}$ be an orthonormal basis for the eigenspace corresponding to the eigenvalue $\lambda_{j}^{2}$. Then, via the Cauchy integral formula, we have

[TABLE]

Similarly, we have

[TABLE]

and

[TABLE]

Combining these terms and noticing that $I_{3}$ cancels with one of the terms in $I_{2}$, we obtain

[TABLE]

To compute the variation of the even wave trace in particular, set $g(z)=\cos(t\sqrt{z})$ in equation (11). Despite the square root, this is in fact an entire function since cosine is even. We have

[TABLE]

Writing this in terms of the wave kernels, we obtain

[TABLE]

Note that all but one of the terms above contain $\dot{\rho}$. Recalling that in the begining of the section, we defined

[TABLE]

to be the coefficient of $\dot{\rho}$ in the expression (12), we obtain Theorem 3.7. Similarly, setting $g(z)=z$ in equation (10) easily yields Theorem 3.8 on the variation of simple eigenvalues. ∎

4. Billiards

Before obtaining a singularity expansion for the wave trace, we first review the relevant background needed on billiards. This will also be useful in our discussion of Chazarain’s parametrix in Section 5.2. In this section, we drop the subscript [math] from our domain in Section 3 and let $\Omega$ denote any bounded strictly convex region in $\mathbb{R}^{2}$ with smooth boundary. This means that the curvature of $\partial\Omega$ is a strictly positive function. The billiard map is defined on the coball bundle of the boundary $B^{*}\partial\Omega=\{(q,\zeta)\in T^{*}\partial\Omega:|\zeta|<1\}$, which can be identified with the inward part of the circle bundle $S_{\partial\Omega}^{*}\mathbb{R}{{}^{2}}$ via the natural orthogonal projection map. We can also identify $B^{*}\partial\Omega$ with $\mathbb{R}/\ell\mathbb{Z}\times[0,\pi]$, where $\ell=|\partial\Omega|$ is the length of the boundary. Define

[TABLE]

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|$. We then define

[TABLE]

where a point $\widehat{(x,\xi)}$ is the reflection of $\xi$ through the cotangent line $T_{x}^{*}\partial\Omega$. In otherwords, $\widehat{(x,\xi)}$ has the same footpoint and cotangential component 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^{*}(\partial\Omega)$. Associated to this map is the billiard flow, or broken bicharacteristic flow, which we denote by $\Phi^{t}$.

The times $t_{\pm}^{j}$ are defined inductively by

[TABLE]

and the maps $\beta^{\pm n}$ are defined via iteration. We also define $T_{\pm}^{n}=\sum_{j=1}^{n}t_{\pm}^{j}$ to be the total time of the flow from $(y,\eta)$ to $\beta^{n}(y,\eta)$.

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 map $\mathcal{C}$ onto the total phase space $B^{*}\partial\Omega$ to obtain a smooth closed curve which is invariant under $\beta$. In the case that the dynamics are integrable, these invariant curves are precisely the Lagrangian tori which folliate the phase space. A point $P$ in $B^{*}\partial\Omega$ is $q$-periodic, $q\geq 2$, if $\beta^{q}(P)=P$. We define the rotation number of a $q$-periodic point $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^{*}\partial\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 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}(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^{*}\partial\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 deeper results and a more thorough treatment of general dynamical billiards, we refer the reader to [Tab05], [Kat05], [Pop94] and [PT11].

4.1. Elliptical billiards

From here on, we let $\Omega$ be an ellipse with horizontal major axis, given by the equation

[TABLE]

The eccentricity of $\Omega$ is defined by $E=a^{2}-b^{2}$.

Birkhoff conjectured that the only strictly convex integrable billiard tables are ellipses. Completely integrable means that there exists a folliation of the phase space by invariant submanifolds. In the context of Hamiltonian systems, this can be shown to be equivalent to the existence of a maximal number of Poisson commuting invariants, called first integrals. The (compact) energy level sets of regular values of these invariants can be shown to be diffeomorphic to tori and the leaves of a maximal such foliation are called Lagrangian tori. The Lagrangian tori are naturally parametrized by so called “angle coordinates” while the transversal directions in phase space are then parametrized by “action coordinates”. It is well known that for each $Z\in(-E,0)\cup(0,b]$, the confocal ellipse ($Z>0$) or hyperbola ($Z<0$) given by

[TABLE]

is also a caustic (see Figure 1). A short proof of this can be found using elementary planar geometry in the appendix of [GM79a]. For elliptical caustics, we follow the notation in [KS18] and [DDCRR17] by setting $\lambda^{2}=Z\geq 0$. In the context of [KS18], integrable is taken to mean that the union of all convex caustics has a non-empty interior in $\mathbb{R}^{2}$. Ellipses are both completely integrable and integrable in the sense of [KS18].

In 1822, Poncelet proved the following remarkable theorem:

Theorem 4.1** ([Pon95b], [Pon95a]).**

Given an ellipse, if a primitive periodic billiard trajectory is tangent to a confocal conic section, then all orbits tangent to that caustic are also periodic, have the same periods, and have the same lengths.

Hence, periodic orbits of a given length come in 1-parameter families. Poncelet’s original theorem, as stated in his 1865 treatises [Pon95b] and [Pon95a], was much more general and concerned inscribing and circumscribing polygons about other confocal conic sections. There are several modern proofs of Poncelet’s theorem. Because of this result, we can project confocal ellipses to invariant curves in the coball bundle $B^{*}\partial\Omega$. We will use $\omega(\mathcal{C})$ for the rotation number of such an invariant curve. For ellipses, the rotation number of a periodic point is always the same as the rotation number of its corresponding invariant curve.

Birkhoff proved in [Bir66] that for any smooth, strictly convex domain and any $\frac{p}{q}\in(0,\frac{1}{2}]$ in lowest terms, there exist at least two geometrically distinct periodic orbits with rotation number $\frac{p}{q}$. A strictly convex billiard table is said to be rationally integrable if for each $q\geq 3$, there exists a caustic consisting of periodic points of rotation number $\frac{1}{q}$. In particular, the rotation number of the invariant curve corresponding to such a caustic must be $\frac{1}{q}$. A third version of Birkhoff’s conjecture is that ellipses are the only rationally integrable strictly convex billiard tables. We have seen that there exist many notions of integrability, yet Birkhoff’s conjecture remains open for all of them. However, as mentioned in Section 2, a local version was proven in [KS18].

In [GM79a], it is shown that periodic points of the billiard map for an ellipse are dense in phase space. Given a length $T\in\text{Lsp}(\Omega)$, the associated fixed point set is denoted $F_{T}=\{(q,\zeta)\in B^{*}\partial\Omega:\Phi^{T}(q,\eta)=(q,\eta)\}$. In [GM79a], the authors construct a special sequence of caustics converging to the boundary such that the associated fixed point submanifold in $B^{*}\partial\Omega$ has exactly two connected components, corresponding to forwards and backwards flow:

Proposition 4.2** ([GM79a]).**

Let $T_{0}=|\partial\Omega|$ be the perimeter of the ellipse $\Omega$. Then in every interval $(T_{0}-\varepsilon,T_{0})$, there exist infinitely many lengths $T\in\text{Lsp}(\Omega)$. For all but a finitely many such $T\in\text{Lsp}(\Omega)$, $F_{T}$ is the union of two invariant curves which are mapped to each other by $(q,\zeta)\to(q,-\zeta)$.

The time reversal map $(q,\zeta)\to(q,-\zeta)$ reverses the direction of a closed geodesic, which of course preserves its length and geometry. We call the lengths $T$ in Proposition 4.2 simple, as all periodic trajectories of length $T$ are tangent to a single caustic. Proposition 4.2 will be crucial in evaluating our singularity expansion in Section 6 and allowing us to differentiate the constants $c_{j}$ and $\widetilde{c_{j}}$ from Theorem 1.4 near $\lambda=0$. That there is only one connected component up to symmetry will rule out any cancellation between terms in the variation of the wave trace.

5. A parametrix for $S_{R}$ and Singularity Expansion

In this section, we use microlocal analysis to obtain a singularity expansion for the variation of the wave trace near the length spectrum. In particular, we microlocalize the wave kernels near periodic transversal reflecting rays in order to obtain a Fourier integral operator (FIO). Using Theorem 3.7, we can rewrite

[TABLE]

where $r_{1},r_{2}$ are the boundary restriction operators, $\Delta:\partial\Omega\to\partial\Omega\times\partial\Omega$ is the diagonal embedding, $\pi_{*}$ is integration over the fibers, and the modified propagator $L_{R}^{b}$ is given by

[TABLE]

This manipulation of notation is just to illustrate how one can decompose the wave trace variation into the composition of simpler Fourier integral operators. However, the wave kernels $L_{R}^{b}$ and $S_{R}$ in our formula are not FIOs near the glancing set $S^{*}\partial\Omega$. In [HZ12], the authors microlocalize the wave kernels near periodic nonglancing orbits and calculate the principal symbol of the composition (13) for Dirichlet and Neumann boundary conditions using the symbol calculus in [DG75]. In contrast to the methods employed in [HZ12], we instead take a more direct approach which avoids an application of the trace formula in [DG75].

We begin by reviewing FIOs and Chazarain’s parametrix for the wave propagator. In Section 5.3, we then cook up an explicit oscillatory integral representation for each term in Chazarain’s parametrix, which microlocally approximates the wave propagator in the interior by using action angle coordinates for the billiard map. This will rely heavily on the symbol calculus in Section 5.2 and a new phase function for wave propagator. This technique can also be extended to deal with other convex billiard tables, although the dynamics are not as simple.

5.1. Fourier Integral Operators

Let $X$ and $Y$ be open sets in $\mathbb{R}^{n_{X}}$ and $\mathbb{R}^{n_{Y}}$ respectively. Recall that a continuous linear operator $A:C_{0}^{\infty}(Y)\to\mathcal{D}^{\prime}(X)$ has an associated Schwartz kernel $K_{A}\in\mathcal{D}^{\prime}(X\times Y)$. 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 form

[TABLE]

is called a Lagrangian or Fourier integral distribution on $X$. If we assume that $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 actually a local diffeomorphism from $C_{\Theta}$ onto $\Lambda_{\Theta}$. In this case, we say that “$\Theta$ parametrizes $\Lambda_{\Theta}$.” The canonical relation or wavefront relation of $A$ is defined by

[TABLE]

which describes how the operator propagates the singularities of distributions on which it acts. More invariantly, one can 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$ given in terms of the parametrization $\iota_{\Theta}$:

[TABLE]

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 Leray measure 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 5.4. For a more thorough reference on the global theory of Lagrangian distributions, see [Dui96]. The order of a Fourier integral operator is defined in such a way that when two Fourier integral operators’ canonical relations meet transversally, then the composition is again a Fourier integral operator and order of the composition is the sum of the orders:

[TABLE]

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 pseudodifferential operators, where $X=Y$ and $m=\mu$ coincides with the order of the corresponding symbol class. A sufficient condition which guarantees that the composition exists is clean or transversal intersections of the two operators’ canonical relations. In general, composition of Fourier integral operators and the associated symbol calculus is somewhat complicated, but is discussed in [Dui96], [Hör71] and [DH72]. We will not directly use the composition formula in what follows.

5.2. Chazarain’s parametrix

Chazarain’s parametrix provides a microlocal description of the wave kernels near periodic transversal reflecting rays. The parametrices for $E_{R}$ and $S_{R}$ are constructed in the ambient Euclidean space $\mathbb{R}\times\mathbb{R}^{n}\times\mathbb{R}^{n}$. We only consider $S_{R}$, as the formula for $E_{R}$ is easily obtained from that of $S_{R}$ by differentiating in $t$. Following the work in [Cha76] and [GM79b], we can find a Lagrangian distribution

[TABLE]

which approximates $S_{R}(t,x,y)$ 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 (15) is locally finite. We first explain what is meant by approximating $S_{R}(t,x,y)$ “microlocally away from the tangential rays.” In general, two distributions $f,g\in\mathcal{D}^{\prime}(\mathbb{R}^{n})$ are said to agree microlocally near a closed cone $\Lambda\subset T^{*}\mathbb{R}^{n}$ if $WF(u-v)\cap\Lambda=\emptyset$. Similarly, using the language from Section 5.1, two operators $A,B:C^{\infty}(Y)\to C^{\infty}(X)$ are said to agree microlocally near a given closed cone $\widetilde{\Lambda}\subset T^{*}X\times T^{*}Y$ if $WF^{\prime}(A-B)\cap\widetilde{\Lambda}=\emptyset$. This second notion is what we will use to say that our parametrix approximates $S_{R}$ microlocally near the canonical relations $\Gamma_{\pm}^{j}$.

In this section, we study the problem

[TABLE]

for $x_{0}\in\text{int}\Omega$ and $K$ a Robin function as in equation (1). Chazarain’s construction begins with the solution of the homogeneous wave equation in the ambient Euclidean space:

[TABLE]

In this case, we have an explicit representation for the fundamental solution of $\Box=\partial_{t}^{2}-\Delta$, given by Kirchhoff’s formula. If we restrict back to $\Omega$, it is clear that the fundamental solution on $\mathbb{R}^{2}$ will in general not satisfy the Robin boundary condition. However, finite speed of propogation implies that it does satisfy the boundary condition for small time, since $u$ vanishes identically in a neighborhood of $\partial\Omega$. If we let $d=d(x_{0},\partial\Omega)$, then we obtain a solution $u$ of the boundary value problem for $|t|<d$. The idea in [Cha76] is to use the billiard flow and properties of FIOs to inductively extend the time interval on which the fundamental solution of the boundary value problem is defined.

With $t_{\pm}^{j}$ defined as in Section 4 on billiards and $(y,\eta)\in T^{*}\Omega$ or $T_{\partial\Omega}^{*}\mathbb{R}^{2}$ inward pointing, we define

[TABLE]

Notice that if $y\in\partial\Omega$ and $|\eta|=1$ is inward pointing, then $\beta(y,\iota^{*}\eta)=\widehat{\lambda_{+}^{1}(y,\eta)}$, where $\beta$ is the billiard map on $B^{*}\partial\Omega$ and $\iota:B\partial\Omega\to S_{\partial\Omega}\mathbb{R}^{2}$ is the natural projection mapping the ball bundle of the boundary to the inward pointing portion of the circle bundle with footpoints on the boundary. Even though we can obtain one from the other, we define both $\lambda_{\pm}$ and $\lambda_{\pm}^{-1}$ in order to separate the forwards and backwards wave propagators corresponding to $t>0$ and $t<0$. After a reflection at the boundary, we can similarly define $\lambda_{\pm}^{j}(y,\eta)$ for any $j\in\mathbb{Z}$. Recall that in Section 4 we defined $T_{\pm}^{j}=\sum_{k=1}^{j}t_{\pm}^{k}$ for $j>0$ and $T_{\pm}^{j}=\sum_{k=j}^{-1}t_{\pm}^{k}$ for $j<0$.

To study how the fundamental solution behaves at the boundary, 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 [Cha76] that the phase functions corresponding to the FIOs $S_{j}$ parametrize the canonical relations

[TABLE]

Again, $j>0$ and $j<0$ correspond to reflections in forward and backward time. In fact, there are four modes of propagation, corresponding to $\pm\tau\geq 0$ and $\pm j\geq 0$ in the canonical relations $\Gamma_{\pm}^{j}$. Chazarain actually showed that there exists FIOs $S_{j}$ such that the sum in (15) is in fact a parametrix for the wave propagator $S_{R}$ with canonical relation

[TABLE]

However, the principal symbols of the operators $S_{j}$ are never computed in [Cha76] and we concern ourselves with the task of explicitly computing them in the special case of an ellipse for the remainder of this section.

Recall that the Hadamard type variational formula for the wave trace in Theorem 3.7 involved the integral of wave kernels over the diagonal of the boundary. To understand the principal symbol first in the interior, we study how the propagator reflects at the boundary. In particular, we want to study the canonical relations $\Gamma_{\pm}^{j}$ restricted to the fibers over the boundary, which we now describe. Denote by

[TABLE]

the fibers corresponding to zero reflections. If we flow out from $A_{\pm}^{0}$ by the Hamiltonian flow $\psi^{\pm t}$ of $H=\tau\pm|\xi|$, we obtain $\Gamma_{\pm}^{0}$. Note that $\psi^{\pm t}$ consists of geodesics lifted to $T^{*}(\mathbb{R}\times\Omega)$. Consider the following subsets of $\Gamma_{\pm}^{0}$:

[TABLE]

If   $\widehat{}$   denotes reflection in the left factor, we have

[TABLE]

We define $A_{\pm}^{j}$ and $\widehat{A_{\pm}^{j}}$ similarly and note that these are precisely the fibers of the canonical relations $\Gamma_{\pm}^{j}$ lying over the boundary. In the next section, we will compute the principal symbol of the wave propagator in coordinates on the critical set. Therefore, we first need to better understand the forwards and backwards symbols on $\Gamma$.

Proposition 5.1**.**

Let $e_{\pm}$ denote the principal symbol of $\widetilde{S_{R}}$ on $\Gamma=\bigcup_{j\in\mathbb{Z},\pm}\Gamma_{\pm}^{j}$. Then, we have

[TABLE]

*Furthermore, the principal symbol for the wave propagator with $K=0$ (Neumann boundary conditions) coincides with $e_{\pm}$ (Robin boundary conditions). *

Proof.

As in [Cha76] and [HZ12], denote by $\sigma_{0}$ the symbol of the restriction to $t=0$, $\sigma_{r}$ the symbol of the boundary restriction operator $r$, and $\sigma_{B}$ the symbol of $rN-K\in\Psi^{1}$. Here, $N$ is an extension of the unit normal vector field on $\partial\Omega$ to a tubular neighborhood of the boundary acting as a differential operator and $K$ is again a Robin function. We have the following implications:

[TABLE]

The first assertion in (16) follows from Theorem 5.3.1 of [DH72] and the remaining formulas are clear. At the boundary, the symbol $\sigma_{B}$ is given by

[TABLE]

on $\Gamma_{\partial\Omega}\circ A_{\pm}^{j}$ and

[TABLE]

on $\Gamma_{\partial\Omega}\circ\widehat{A}_{\pm}^{j}$. The symbol of $rK$ doesnt appear since multiplication by $K$ is a $\Psi$DO of order [math] while, $N\in\Psi^{1}$. Hence, the fourth equation in (16) implies that on the boundary, we have

[TABLE]

Note that in equation (17), both sides involve the composition $\sigma_{r}\circ e_{\pm}$. The formula for the principal symbol of the composition of FIOs is quite complicated, but is discussed more thoroughly in [GM79b], [HZ12], [DG75] and [DH72]. Since $\langle\lambda(y,\eta),\nu_{y}\rangle$ is nonvanishing, equation (17) tells us that the direct and reflected symbols coincide on the boundary. Multiplying the second equation in (16) by $\tau$ and adding/subtracting it to the third equation gives

[TABLE]

in the interior. It is elementary to see that $\sigma_{Id}$ is the canonical half density $|dt\wedge dy\wedge d\eta|^{1/2}$. The first equation in (16) implies that the symbol is invariant under geodesic flow, so the claim follows on $\Gamma_{\pm}^{1}$. In fact, since we have already noted that the direct and reflected symbols coincide over the boundary, the fact that $\Gamma_{\pm}^{j}$ is the flowout of $\widehat{A}_{\pm}^{j}$ then implies that the claim extends to all $\Gamma_{\pm}^{j}$. ∎

Remark 5.2*.*

Chazarain’s parametrix actually computes the full symbol by solving successive transport equations and Borel summing the terms. The full proof and explicit computation of the symbol can be found in the original French paper [Cha76]. The actual solution operator $S_{R}$ could be obtained from $\widetilde{S_{R}}$ by adding correction terms via Duhamel’s principle. However, we only need the principal symbol in our calculation. In [GM79b], a more general situation is treated in which both $M=\Omega\times\{t=0\}$ and $M^{\prime}=\partial\Omega\times\mathbb{R}$ are nonglancing, noncharacteristic hypersurfaces for the wave propagator. A Fourier integral operator is then constructed iteratively to solve the localized hyperbolic pseudodifferential equation.

We now make precise the notion of microlocalized FIOs. Recall Theorem 4.2 in Section 4, which provides an ample number of caustics having simple length in any neighborhood of $|\partial\Omega|$. Hence, for $j\in\mathbb{Z}$ large and positive, we can consider periodic orbits having simple length $T_{j}$, making a single rotation and precisely $j$ reflections at the boundary. We would like to microlocalize $S_{R}$ near orbits of such a simple length $T_{j}$. Let $\chi_{1}(t)$ be a smooth cutoff function which is identically equal to $1$ on an open neighborhood of $T_{j}$ and vanishes in a neighborhood of all other $T\in\text{Lsp}(\Omega)$. As we remarked above, 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 $\cup_{\pm}\Gamma_{\pm}^{j}$ and is conic in the fiber variables $\tau,\xi$ and $\eta$. Quantizing $\chi_{2}$ gives a pseudodifferential operator with wavefront set contained in support of $\chi_{2}$. For a reference, see Chapter 18 of [Hör85a]. 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_{R}$ is then smoothing away from the periodic orbits of rotation number $1/j$. Since $T_{j}$ was assumed to be simple, the trace of the above composition is equal to the wave trace modulo $C^{\infty}$ in a neighborhood of $T_{j}$.

5.3. Computing the singularity in elliptical polar coordinates

In the previous section, we reviewed Chazarain’s parametrix and computed the principal symbol and canonical relation for the wave propagator. In contrast to the methods employed in [HZ12], we now want to cook up an oscillatory integral such that microlocally near $\Gamma_{\pm}^{j}$,

[TABLE]

where $S_{j}$ is the $j^{th}$ term in Chazarain’s parametrix corresponding to a wave with $j$ reflections. Here, $L.O.T$ denotes lower order terms in the sense of Lagrangian distributions. Due to the presence of different Maslov factors for $\pm\tau>0$ (see Sections 5.1 and 5.4), it is actually more convenient to find operators

[TABLE]

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 find suitable phase functions $\Theta$ parametrizing $\Gamma_{\pm}^{j}$, which we can do only after learning more about elliptical billiards. The following geometric description of almost periodic orbits in the ellipse is crucial:

Lemma 5.3**.**

For $j\in\mathbb{Z}$ sufficiently large and any two points $x,y\in\text{int}\Omega$ near 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.

We first explain what is meant by approximately one roation. Let $\xi\in S_{x}^{*}\Omega$ be one of the $4$ covectors corresponding to the initial condition of a counterclockwise orbit described in Lemma 5.3. 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 before 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 7). 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 Lemma 5.3, we have

[TABLE]

Here, $\ell=|\partial\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 4. The choice of $\ell/100$ is somewhat arbitrary, but having $\ell$ in the numerator allows for scale invariance and finding the optimal constant in the denominator is irrelevant for our purposes. The notion of approximately one clockwise rotation is defined similarly. We relegate the proof this theorem to Section 7, as it uses formulas we haven’t yet discussed and is of independent interest. The proof of Lemma 5.3 actually provides more information. Of the four counterclockwise orbits emanating from $x$, two of them become tangent to a confocal ellipse 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 confocal ellipse and we call these $N$ orbits (for nontangency) with first link 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 confocal ellipse (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 counterclockwise orbits from $x$ to $y$, which we denote by $TT$, $TN$, $NT$, and $NN$. See Figure 2 for an example with $j=4$. The same characterization also applies to the four clockwise orbits in Lemma 5.3, which can be obtained by reflecting the domain through the vertical axis, finding all orbits making approximately one counterclockwize rotation from $x$ to $y$, and then reflecting these orbits back through the vertical axis (see Section 7 for a more detailed discussion). These configurations will be important in determining which limiting orbits give periodic trajectories of precisely $j$ reflections as $(x,y)\to\Delta\partial\Omega$.

Definition 5.4**.**

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 Lemma 5.3. It depends only on $x,y,j$ and $k$. We use the convention that the indices $1\leq k\leq 4$ correspond to the counterclockwise orbits $TT,TN,NT,NN$ and the indices $5\leq k\leq 8$ correspond to their reflections about their clockwise counterparts (reflections of the first four orbits through the vertical axis).

The author learned of a similar function in [MM82] (page 492), where its restriction to the boundary is defined. In such a case, i.e. if $x,y\in\partial\Omega$, it is shown in [GM81], [MM82] and [Pop94] that only a single counterclockwise orbit of $j$ reflections exists between the boundary points if they are sufficiently close and $j$ is sufficiently large. Upon inspection of the geometric proof given in Section 7, one can actually see that as $x$ and $y$ approach the diagonal of the boundary from the interior, the corresponding orbits coalesce and converge to the orbits described in [MM82]. However, the limiting orbits may have a different number of reflections (see proof of Lemma 5.8). We define phase functions $\Theta_{j,\pm}^{k}$ by the formula

[TABLE]

Lemma 5.5**.**

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 graphs $\Gamma_{\pm}^{j}$. In particular, both $\Gamma_{+}^{j}$ and $\Gamma_{-}^{j}$ are unions of $8$ connected components, which we denote by $\Gamma_{\pm}^{j,k}$.

Proof.

For any $x,y\in\Omega$ let

[TABLE]

denote the length functional. We first show that billiard trajectories from $x$ to $y$ are in one to one correspondence with critical points of (18) with respect to $q\in\partial\Omega^{j}$. Let $g\in C^{\infty}(\mathbb{R}^{2})$ be a defining function for $\partial\Omega$ and consider $q$ as a variable in $\mathbb{R}^{2}\times\cdots\times\mathbb{R}^{2}=\mathbb{R}^{2j}$ rather than $\partial\Omega^{j}$. If $q$ is a critical point of (18), then as in the method of Lagrange multipliers, by setting $x=q_{0}$ and $y=q_{j+1}$, we find that for $1\leq m\leq j$, there exists $\lambda_{m}\in\mathbb{R}$ such that

[TABLE]

Since $\nabla_{q_{m}}g\perp\partial\Omega$, this implies that the two unit vectors in the formula for $\partial_{q_{m}}L_{x,y}$ have opposite tangential components, which is precisely the condition giving elastic collision at the boundary (angle of incedince equals angle of reflection). Similarly, if this condition is satisfied, then $q$ is a critical point for (18).

We now consider the functions $\Psi_{j}^{k}$ in Definition 5.4. We have

[TABLE]

where $q_{m}^{k}(x,y)$ is the $m$th impact point on the boundary for the billiard trajectory corresponding to $\Psi_{j}^{k}$. As opposed to the $q_{m}$ in the length functional (18), $q_{m}^{k}$ will in general have a nontrivial dependence on $x$ and $y$. Differentiating (19) in $x$, we obtain

[TABLE]

Since for each $x,y\in\Omega$, the path defined by $(x,q^{k},y)$ corresponds to a billiard trajectory, we see that all of the terms except the first telescope in (20). Hence,

[TABLE]

Similarly, differentiating (19) in $y$, we obtain

[TABLE]

Geometrically, these gradients are the incident and (reflected) outgoing unit directions of the billiard trajectories described in Lemma 5.3.

We now consider the maps

[TABLE]

on the critical set $C_{\Theta_{j,\pm}^{k}}=\{t-\Psi_{j}^{k}=0\}$. Inserting formulas (21) and (22) into (23) and comparing with the canonical graphs

[TABLE]

from Section 5.2, we see that $\iota_{{\Theta_{j,\pm}^{k}}}:C_{\Theta_{j,\pm}^{k}}\to\Gamma_{\pm}^{j}$ is a local diffeomorphism. Since $1\leq k\leq 8$, it follows that both $\Gamma_{+}^{j}$ and $\Gamma_{-}^{j}$ are the unions of $8$ connected components. ∎

We now want to derive an explicit formula for the principal symbol $e_{\pm}$ of $S_{R}$ in coordinates. Referring to formula (14) for the principal symbol of an FIO, we easily see that in our setting, $dC_{\Theta_{j,\pm}^{k}}=\mp d\tau\wedge dx\wedge dy$. In Proposition 5.1, we calculated that $e_{\pm}=\frac{1}{2\tau i}|dt\wedge dy\wedge d\eta|^{1/2}$. Since we now know that the phase functions $\Theta_{{j,\pm},\pm}^{k}(t,\tau,x,y)=\pm\tau(t-\Psi_{j}^{k}(x,y))$ locally parametrize the connected components of $\Gamma_{\pm}^{j}$, we now want to calculate $a_{0}$ (as in Section 5.1) by changing variables. It is ultimately more convenient to introduce a conformal change of coordinates which is suitable to computing the symbol $e_{\pm}$ for the ellipse:

Definition 5.6**.**

Elliptical polar coordinates are defined on $\mathbb{R}^{2}\backslash[-c,c]$ by the equations:

[TABLE]

Here, $(\mu,\varphi)\in(0,\infty)\times\mathbb{R}/(2\pi\mathbb{Z})$ and $c=\sqrt{a^{2}-b^{2}}$ is the semifocal distance, i.e. the distance between the origin and a focal point of the ellipse

[TABLE]

It is easy to check that the Euclidean metric in these coordinates is a conformal multiple of $d\mu d\varphi$:

[TABLE]

In particular, the vector fields $\partial/\partial\mu$ and $\partial/\partial\varphi$ are orthogonal at each point. For a fixed $\mu>0$, the $\varphi$ coordinate parametrizes a confocal ellipse of eccentricity $1/\cosh(\mu)$. When projected onto the phase space $B^{*}\partial\Omega$, these curves are precisely the invariant Lagrangian tori for the billiard map. In particular, fixing

[TABLE]

gives a parametrization of the boundary $\partial\Omega$. Similarly, for a fixed $\varphi>0$, the $\mu$ coordinate parametrizes a branch of a confocal hyperbola. The reason we use these coordinates is because up to a conformal factor, tangential differentiation on the ellipse becomes $\partial/\partial\varphi$ while differentiating in the normal direction becomes $\partial/\partial\mu$. Since the wave kernel $S_{R}$ is a function of both the $x$ and $y$ variables, we use elliptic coordinates for $y$ as well:

[TABLE]

Note that the Leray form is coordinate independent. Hence, without loss of generality, we compute that in elliptical coordinates,

[TABLE]

On the critical set, we have

[TABLE]

From now on, we drop the the $j,k$ subscripts and write $\Psi_{j}^{k}=\Psi$ so that we may use subscripts to denote derivatives. We have

[TABLE]

Wedging all these terms, we find that

[TABLE]

Keeping in mind that all $\Psi$ terms depend on $j$ and $k$, we denote the above factor by

[TABLE]

Then, on each of the canonical relations $\Gamma_{\pm}^{j,k}$ (see Lemma 5.5), we have

[TABLE]

As a result, we have proved:

Theorem 5.7**.**

Microlocally near $\Gamma_{\pm}^{j}$, the following oscillatory integral is a parametrix for $S_{j}$ in an open neighborhood of $\Delta\partial\Omega\subset\Omega\times\Omega$:

[TABLE]

We denote the operators in this sum by $S_{j,\pm}^{k}$ and and also define $S_{j}^{k}=S_{j,+}^{k}+S_{j,-}^{k}$.

Recall that according to Theorem 3.7, the variation of the wave trace is given by

[TABLE]

where $L_{R}^{b}$ is defined as

[TABLE]

The highest order terms come only from the differentiated sine kernels, which implies that

[TABLE]

In fact, we can discard even more terms in the singularity expansion:

Lemma 5.8**.**

Modulo Maslov factors and distributions of lower order, the variation of the localized (even) wave trace near a simple length $T_{j}$ is given by

[TABLE]

Proof.

For the localized wave trace, we only need to consider orbits which contribute to the singularity at $T_{j}$. Recall that for positive time, Lemma 5.3 gives $8$ orbits connecting $x$ to $y$. These orbits coalesce into one of the orbits from [MM82] as $(x,y)\to\Delta\partial\Omega$. However, as the orbits coalesce within various configurations, not all of the limiting orbits will have $j$ reflections. As $T_{j}$ is simple, only the limiting periodic orbits having exactly $j$ reflections will contribute to the wave trace near $t=T_{j}$. Figure 2 may be useful in visualizing the geometric arguments which follow. As $(x,y)\to\Delta\partial\Omega$, the two corresponding orbits in $TT$ configuration $(k=1,5)$ converge geometrically to a periodic orbit of $j+1$ reflections. The additional vertex appears at the boundary point where $x$ and $y$ coalesce. Similarly, the $NN$ orbits $(k=4,8)$ can be seen to converge to a periodic orbit 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 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 the boundary, only $4$ out of the $8$ orbits contribute to periodic trajectories of $j$ reflections on the boundary. However, in the limit, two additional $TT$ orbits of $j-1$ reflections converge to a periodic orbit of $j-1+1=j$ reflections. Similarly, two $NN$ orbits of $j+1$ reflections converge to a periodic trajectory 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$ impact points at the boundary cannot converge to a periodic orbit of $j$ reflections. As we have localized the wave trace near the simple length $T_{j}$, only the $4+2+2=8$ orbits which converge geometrically to a periodic orbit of exactly $j$ reflections will contribute to the singularity near $T_{j}$. All additional orbits contribute smooth errors to the wave trace in a small neighborhood of $T_{j}$. It should also be clarified that although the parametrices $S_{j}$ 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_{R}(t,q,q)$ 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}$ in fact shows that they extend continuously up to the boundary since the functions $\Psi_{j}^{k}$ do. The true wave kernels $S_{R}$ also extend continuously up to the boundary as a family of distributions. To see this, note that

[TABLE]

where $(\psi_{j})_{j=1}^{\infty}$ is an $L^{2}$ orthonormal basis of Robin eignenfunctions. Multiplying by a test function $\varphi$ in time and integrating by parts $4k$ times, we see that

[TABLE]

Combining Weyl’s law on the asymptotic growth of $\lambda_{j}$ (see [Zay04], [Ivr16]) and Hörmander’s $L^{\infty}$ bounds for eigenfunctions ([Hör68]), we see that the integrand in (25) can be made absolutely convergent for $k$ sufficiently large. An application of the dominated convergence theorem then shows that (25) is actually smooth in $x,y$, so $S_{R}$ 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, which is required for the trace formula. ∎

Definition 5.9**.**

As shown in the proof of Lemma 5.8 above, for each $j$, there exist $8$ limiting trajectories which converge geometrically to periodic orbits 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}$. By the results in [MM82], [GM79a] and [Pop94], 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}$ actually coincide for $x,y\in\partial\Omega$ near the diagonal. We denote their common value by $\Psi_{j}$.

As we obtained a rather explicit formula for $S_{j}$ in Theorem 5.7, it now remains to differentiate the kernels $S_{j}^{k}$ and substitue them into Lemma 5.8. Using our oscillatory integral representation for $S_{j,\pm}^{k}$ in Theorem 5.7, we find that microlocally near $\Gamma_{\pm}^{j,k}$ and $t=T_{j}$,

[TABLE]

in an open neighborhood of the diagonal of the boundary. We have only written the terms coming from $-\nabla^{T}\nabla^{T}-\Delta$ acting on the phase function in equation (26), as all other terms don’t contribute positive powers of $\tau$ and can be regarded as lower order in the singularity expansion. The operator $\nabla^{\perp}$ in the integrand of (26) is a conformal multiple of the vector field $\frac{\partial}{\partial\nu}$ coming from elliptical polar coordinates, which gives an extension of the normal vector field to a neighborhood of the boundary. As Theorem 3.7 tells us that the variation of the wave trace is given by integrating the kernels $L_{R}^{b}$ and $S_{R}$ over the diagonal of the boundary, we want to understand the restriction of (26) to the boundary. In Definition 5.9, we noted that the length of the unique orbit connecting two boundary points with $j$ reflections is well defined. When $x=y=q\in\partial\Omega$, $\Psi_{j}(q,q)$ gives precisely the length $T_{j}$ of a periodic orbit with $j$ reflections emanating from and terminating at $q$. By Poncelet’s Theorem (4.1), $\Psi_{j}(q,q)$ is actually equal to the constant function $T_{j}$, which simplifies the phase in equation (26). The differentiated kernels in equation (26) also have factors of $A_{j}^{k}$ and $\nabla_{2}^{\perp}\Psi_{j}^{k}$ in the integrand. We now discuss how to extend these derivatives of $\Psi_{j}^{k}$ to the boundary in a manner analogous to that of Definition 5.9.

We have already computed the $x,y$ gradient of the functions $\Psi_{j}^{k}$ in equations (21) and (22) in the proof of Lemma 5.5:

[TABLE]

Geometrically, these are the incident and reflected outgoing unit directions of the corresponding billiard trajectories at $x$ and $y$. The expression $|\nabla_{y}^{\perp}\Psi_{j}^{k}|^{2}$ in (26) can easily be seen to be $\sin^{2}\omega_{j}^{k}$, where $\omega_{j}^{k}$ is the angle made between the terminal link of the billiard trajectory and the positively oriented tangent line to the confocal ellipse on which $y$ lies. As $x,y\to\partial\Omega$, the absolute value of these angles associated to trajectories in the $\mathcal{G}_{j}$ converge to the terminal angle of the unique limiting orbit connecting boundary points in [MM82]. We are careful to point out that only the absolute values of the angles converge, since the angles associated to orbits in $TN$ and $NN$ configurations actually converge to minus the angle of incidence of the limiting trajectory. All of the limiting orbits which connect boundary points in [MM82] are automatically in $TT$ configuration.

In order to understand the factor $A_{j}^{k}$, we must compute the Hessian of $\Psi_{j}^{k}$ and its restriction to the boundary. Recall that $A_{j}^{k}$ is given by equation (24) in elliptical polar coordinates. We can further simplify that expression to obtain the following:

Lemma 5.10**.**

On the boundary, all of 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}^{8}$ corresponding to orbits in $\mathcal{G}_{j}$ coincide up to sign. We denote their common (absolute) value by $|A_{j}|$, which on $\Delta\partial\Omega$ in particular, satisfies the equation

[TABLE]

Here, $\omega$ is the angle of incidence of the unique periodic orbit with $j$ reflections at $(a\cos\theta,b\sin\theta)\in\partial\Omega$ and $f(\mu_{0},\theta)={({a^{2}(\cosh^{2}\mu_{0}-\cos^{2}\theta)})}^{1/2}$ is the inverse of the conformal factor in elliptical coordinates.

Proof.

Let $(m,k)$ denote an admissible pair of indices corresponding to an orbit $\gamma_{m,k}\in\mathcal{G}_{j}$. Recall the notation in the proof of Lemma 5.5, where we described a billiard trajectory by the point $(x,q,y)\in\Omega\times\partial\Omega^{m}\times\Omega$. If $x\in\Omega$, let us denote the angle between $x-q_{1}$ and the positively oriented tangent line to the confocal ellipse on which $x$ lies by $\omega_{1}$. Similarly, if $y\in\Omega$, let us also denote the angle between $y-q_{j}$ and positively oriented tangent line to the confocal ellipse on which $y$ lies by $\omega_{2}$. Since the elliptical coordinates $(\mu,\varphi,\nu,\theta)$ are conformally equivalent to Euclidean $(x,y)$ coordinates, we have

[TABLE]

Equations (21) and (22) in the proof of Lemma 5.5 then tell us that

[TABLE]

where the $\pm$ in the equations for $\Psi_{\mu}$ and $\Psi_{\nu}$ are dependent on the configuration of the orbit. They are $+$ if the corresponding initial or final link is a $T$ link and $-$ if it is an $N$ link. Using (27) to calculate the second derivatives, we have

[TABLE]

Then, inserting (28) into the expression (24) for all possible configurations, we find that on the boundary,

[TABLE]

where the superscripts on $\mathcal{G}_{j}$ indicate the subcollection of orbits within a particular configuration. Before evaluating this expression on the diagonal of the boundary, we differentiate $\omega_{1}$ in the direction $L$ of the last link to see that

[TABLE]

where the $\pm$ correspond to whether the last link is a $T$ link $(+)$ or an $N$ link $(-)$. This implies that,

[TABLE]

Note that on the diagonal of the boundary, $\omega_{1}=\omega_{2}$ by the law of equal reflection for billiards. We denote their common value by $\omega$. Inserting formula (30) into (29) and evaluating on $\Delta\partial\Omega$, we find that

[TABLE]

which proves the lemma. ∎

Substituting the formula for $A_{j}$ in Lemma 5.10 into (26) and performing the integral (26) in $\tau$, we obtain:

Corollary 5.11**.**

The variation of the wave trace localized near $t=T_{j}$ is

[TABLE]

where $\sigma$ is a Maslov index.

Proof.

For a fixed $j$, the parametrices corresponding to $\Gamma_{+}^{j}$ and $\Gamma_{-}^{j}$ are first multiplied by Maslov factors of the form $e^{\pm i\sigma\pi/4}$. This is due to the multiplicity of phase functions parametrizing the canonical relation of the wave propagator, as briefly described in Section 5.1. It is well known that the Maslov factors on the two branches $\Gamma_{-}^{j}$ and $\Gamma_{+}^{j}$ of the canoncial relation are conjugate to one another, owing to the two modes of propagation:

[TABLE]

While in principle, the Maslov indices might depend on $j$, it is shown in [GM79b] that the Maslov indices in Chazarain’s parametrix remain unchanged after a reflection at the boundary. The Maslov indices $\pm\sigma$ will in fact be explicitly computed in Section 5.4 below. The contributions of the wave kernels corresponding to $\Gamma_{+}^{j}$ and $\Gamma_{-}^{j}$ are then added and the result follows from a limiting argument for the Fourier transform of the homogeneous distribution $\tau^{3/2}\mathbbm{1}_{(0,\infty)}$, as can be found in Chapter 7 of [Hör03]. Since the $8$ terms in (26) corresponding to orbits in $\mathcal{G}_{j}$ coincide on the boundary, we multiply the final integral by a factor of $8$. ∎

5.4. Computing the Maslov index

To explicitly compute the Maslov factors $e^{i\pi\sigma_{j}^{\pm}/4}$ on $\Gamma_{\pm}^{j}$, we use an argument due to Keller ([Kel58]), following the presentation in [FHH18a]. The free wave propagator $U(t)=e^{-it\sqrt{-\Delta}}$ on $\mathbb{R}^{2}$ has an integral kernel given by

[TABLE]

considered as a distributional half density (see Section 5.1). By changing variables and applying the method of stationary phase, it is shown in [FHH18a] that the principal symbol of $U(t)$ on $N^{*}\{|x-y|=t\}=\Gamma_{\pm}^{0}$ is

[TABLE]

Hence, the Maslov indices are given by $\sigma_{0}^{\pm}=\pm 1$ on $\Gamma_{\pm}^{0}$. As mentioned in the proof of Corollary 5.11, the arguments in [GM79b], which are in turn based on the construction in [Cha76], 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}$. Since both the kernels on $\Gamma_{\pm}^{j}$ contribute to the wave trace singularity near $T_{j}$, we sum together the contributions of $S_{j}^{-}$ and $S_{j}^{+}$ which explains the real parts in Theorem 1.4 and Corollary 5.11.

5.5. Calculating $\frac{\partial\omega}{\partial\theta}$

To calculate the angular derivative in Lemma 5.10, we will first relate it to the billiard map and then utilize some special dynamical properties of the ellipse in action angle coordinates. Recall that the billiard map takes place on the coball bundle $B^{*}\partial\Omega$, which is diffeomorphic to the inward facing portion of the circle bundle with footpoints on the boundary and can be parametrized by coordinates $(\varphi,\omega)\in S^{1}\times S^{1}$ (although we only consider the nontangential, inward pointing directions corresponding to $0<\omega<\pi$). Between any two points $x,y\in\partial\Omega$, the results of [MM82] show that there exists a unique broken geodesic of $j$ reflections emanating from $x$ and terminating at $y$. This geodesic makes an initial angle of $\omega$ (depending on both $x$ and $y$) with the tangent line $T_{x}\partial\Omega$. Setting $x=y$ above gives the angle $\omega(x)$ corresponding to a periodic orbit, which we considered in Lemma 5.10. Letting $\varphi,\theta\in S^{1}$ be the angular variables which parametrize $x$ and $y$ respectively in elliptical polar coordinates, we need to calculate the quantity $\frac{\partial\omega}{\partial\theta}$, evaluated on the diagonal $\{\varphi=\theta\}$. Consider the map

[TABLE]

where $\pi_{1}$ is the projection onto the first factor. Fixing $\varphi$ and differentiating both sides in $\theta$ gives

[TABLE]

Hence, we have

[TABLE]

We will use formula (34) to calculate $\frac{\partial\omega}{\partial\theta}$. Recall that the linearized Poincaré map of the iterated billiard map $\beta^{j+1}$ at a periodic point $(\varphi,\omega)$ is given by

[TABLE]

where $\beta_{1}^{j+1}$ and $\beta_{2}^{j+1}$ are the first and second components (in elliptical coordinates) of $\beta^{j+1}$. We are precisely interested in the $(1,2)$ entry of this matrix.

To evaluate this quantity, we will use action angle coordinates for elliptical billiards, which we now describe, following the presentation in [KS18], [DDCRR17] and [CF88]. We begin by developing some basic elliptic function theory.

Elliptic functions and elliptic integrals were first studied in the context of computing the arclength of an ellipse. It is therefore no surprise that these same objects appear naturally in the study of elliptical billiards. Formally, an elliptic function is given by a doubly periodic, meromorphic function on the complex plane. One way to obtain elliptic functions is by inverting elliptic integrals:

Definition 5.12**.**

An incomplete elliptic integral of the first kind is an integral of the form

[TABLE]

The quantity $\varphi$ is referred to as the amplitude and $k$ is the modulus. A complete elliptic integral of the first kind is given by fixing $\varphi=\pi/2$:

[TABLE]

Note that for a fixed $k$, $F(\varphi,k)$ is an increasing function of $\varphi$. The amplitude function $\text{am}(s;k)$ is obtained by inverting $F$ in the variable $\varphi$:

[TABLE]

Definition 5.13**.**

The Jacobi elliptic functions are defined by

[TABLE]

These are elliptic functions with periods $4K(k)$ and $4iK(k^{\prime})$, where $k\in(0,1)$ is called the modulus and $k^{\prime}=\sqrt{1-k^{2}}$ is the complimentary modulus. The reason these elliptic functions are useful is that they provide coordinates on phase space in which the billiard map becomes a simple translation. To the confocal ellipse

[TABLE]

we associate the following parameters:

[TABLE]

Let us also denote the $4K(k_{\lambda})$ periodic boundary parametrization associated to the caustic $C_{\lambda}$ by $q_{\lambda}(s)=(-a\text{sn}(s;k_{\lambda}),b\text{cn}(s;k_{\lambda}))$ (a reflection about the $y$ axis of the parametrization considered in [DDCRR17]). It is proven in [CF88] that for all $s\in\mathbb{R}$ the line segment connecting $q_{\lambda}(s)$ and $q_{\lambda}(s+\delta_{\lambda})$ is tangent to the caustic $C_{\lambda}$. In other words, $(s,\lambda)$ are precisely the action-angle coordinates from Section 4 and in these coordinates, the billiard map is given by a linear rotation along the invariant tori, which are the projections of $C_{\lambda}$ onto $B^{*}\partial\Omega$.

To relate these action angle coordinates to the elliptical polar coordinates in which we calculated $A_{j}$ (see Lemma 5.10), note that

[TABLE]

which implies

[TABLE]

Here, $t_{\varphi}$ is defined implicitly by the equations (36) and (38). Similarly, we find that

[TABLE]

Differentiating $B_{j}$ in $\theta$, we obtain

[TABLE]

By the chain rule, we see that

[TABLE]

We can factor out $\frac{\partial\lambda^{2}}{\partial\omega}$ from (39) and calculate each of the individual terms explicitly. Using the implicit function theorem, we find that

[TABLE]

evaluated at the point $s=t_{\varphi}+(j+1)\delta_{\lambda}$. We now find the derivative of $t_{\varphi}$ with respect to $\lambda$. Since $x$ is fixed, we know that the argument $\varphi=\text{am}(t_{\varphi};k_{\lambda})$ must also be fixed. Differentiating in $\lambda$ under the integral, we see that

[TABLE]

Using the formula (37) for $\delta_{\lambda}$, we also calculate that

[TABLE]

For the last term in (39), write

[TABLE]

and differentiate both sides in the variable $k_{\lambda}^{2}$. We find that

[TABLE]

where $s$ is evaluated at $t_{x}+(j+1)\delta_{\lambda}$. It is also easy to see that

[TABLE]

At the critical $\lambda$ corresponding to the angle ${\omega}(x)$ generating a periodic orbit at $x$, we have $\text{am}(t_{\varphi}+(j+1)\delta_{\lambda})=\text{am}(t_{\varphi})+2\pi$. Using this, we calculate that

[TABLE]

We have not simplified the expression in parentheses in (41), since what will ultimately be important is that this term depends only on $\lambda$ and not on $\varphi$ or $\omega$. Let us denote this factor by

[TABLE]

The term $\partial\lambda^{2}/\partial\omega$ is more difficult to compute and will rely on a geometric lemma, which we now present. As the ellipse is folliated by caustics, for each $(x,\omega)\in S_{\partial\Omega}^{*}(\mathbb{R}^{2})$, there exists a unique $\lambda$ such that the line segments of the billiard flow are always tangent to the confocal ellipse/hyperbola of parameter $\lambda$. The following lemma expresses this relationship between the confocal caustic and angle of incidence:

Lemma 5.14**.**

The billiard ray emanating from $(a\cos\varphi,b\sin\varphi)$ at angle $\omega$ is tangent to the elliptical caustic $C_{\lambda}$, where the relationship between $\lambda$ and $\omega$ is given by

[TABLE]

Proof.

For simplicity, consider the complexified parametrization of the ellipse given by $\gamma(\varphi)=a\cos\varphi+ib\sin\varphi$. The tangent line at $\gamma(\varphi)$ is then parametrized by

[TABLE]

Rotating this line counterclockwise by the angle $\omega$, we see that the billiard ray is parameteriszed by

[TABLE]

Taking real and imaginary parts, we find that

[TABLE]

For a given $\varphi,\omega$, there exist infinitely many caustics which intersect the line $L_{\omega}$. However, only one such caustic intersects $L_{\omega}$ at a single point of tangency. To find the parameter of this caustic, we look for a solution of the equation

[TABLE]

This is a quadratic equation in the variable $t$ and a caustic corresponding to a point of tangency will give rise to a repeated root. Thus, to find $\lambda^{2}$, we set the discriminant of this equation equal to zero. For convenience let us put

[TABLE]

If we set the discriminant equal to zero, we obtain

[TABLE]

After some obvious cancellations and simplifications, multiplying both sides by $(a^{2}-\lambda^{2})(b^{2}-\lambda^{2})$ gives

[TABLE]

Two simple computations show that

[TABLE]

and

[TABLE]

Plugging these into the above equation for $\lambda^{2}$ completes the proof of the lemma.

∎

Differentiating the formula in Lemma 5.14 by $\omega$ gives

[TABLE]

which is the last term in (39) we needed to calculate.

Remark 5.15*.*

In [GM79b], a different formula is derived in the usual circular polar coordinates with angular parameter $\alpha$. In this case,

[TABLE]

where

[TABLE]

The relationship between $\alpha$ and the elliptical angular coordinate $\varphi$ is given by

[TABLE]

6. Converting the singularity into an elliptic integral

6.1. Proof of $\dot{\rho}=0$

It now remains to convert the singularity expansion in Corollary 5.11 into an elliptic integral, following the ideas in [GM79a]. Recall that according to Corollary 5.11, the variation of the localized wave trace is given by

[TABLE]

The computations in Section 5.5 combined with Lemma 5.10 lead us to the formula

[TABLE]

on the diagonal of the boundary, where $G(\lambda)$ is given by (42) and $\frac{\partial\lambda^{2}}{\partial\omega}$ is given by (43). The formulas (42) and (40) show that $G$ is in fact a nonvanishing analytic function of $\lambda$. Putting this all together, we see that the principal term in the variation of the wave trace in Corollary 5.11 is given by the product of the distribution

[TABLE]

and the factor

[TABLE]

Since we are evaluating on the diagonal of the boundary, Poncelet’s Theorem (4.1) guarantees that the parameter $\lambda$ corresponding to a $j$-periodic geodesic is in fact independent of $\varphi$. Hence, the $G(\lambda)$ factor can be pulled outside of the integral in (45). However, both $\omega$, and $t_{\varphi}$ in the integrand depend on $\varphi$, so it remains to compute this dependency and parametrize the boundary explicitly. Recall that Lemma 5.14 gives us

[TABLE]

Since $\omega\notin\{0,\pi/2,\pi\}$ for orbits which are tangent to a confocal ellipse, this equation determines the trigonometric terms appearing in the integrand of (45) up to a sign:

[TABLE]

To simplify notation, set $C(\varphi)=\left(b^{2}+(a^{2}-b^{2})\sin^{2}(\varphi)\right)$. Recalling that in elliptic coordinates, if we fix

[TABLE]

then we obtain a parametrization of $\partial\Omega_{0}$ in terms of $\varphi\in\mathbb{R}/2\pi\mathbb{Z}$, given by

[TABLE]

In these coordinates, the line element on $\partial\Omega_{0}$ is

[TABLE]

Recalling that $\cos\varphi=-\text{sn}(t_{\varphi};k_{\lambda})$ due to our convention (38) on the phase shift, a simple computation using formula (46) then shows that

[TABLE]

which simplifies the integrand in (45) to

[TABLE]

Reinserting this into the integrand in (45), we obtain

[TABLE]

from which the first formula in Theorem 1.4 follows.

We now set the integral in equation (48) equal to zero. Since $\lambda$ is nonzero and independent of $\varphi$, we can divide out the separated terms depending only on $\lambda$ from both sides and what is left will be a nonzero analytic function $F(\varphi)$ multiplied by $\dot{\rho}(\varphi)(C(\varphi)-\lambda^{2})^{-1/2}$:

[TABLE]

Here, we have written

[TABLE]

Using the $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetry condition on $\dot{\rho}$, the equation (49) reduces to

[TABLE]

We now recall that since the singularity expansion is localized near a simple length $T_{j}$, we have $\lambda=\lambda_{j}$ and $\lambda_{j}\to 0$ as $j\to\infty$. The expression (50) is actually analytic in the parameter $\lambda^{2}$ and vanishes at each $\lambda_{j}^{2}$. Since it has an accumulation point of zeros, it is actually flat at $\lambda=0$. Differentiating $k$ times under the integral (50) in the parameter $\lambda^{2}$ and evaluating at $\lambda=0$, we see that

[TABLE]

It is clear that the functions $|C(\varphi)|^{-k}$ form a subalgebra of $C(S^{1})$. Since we have restricted the domain to $(0,\pi/2)$, this subalgebra also separates points, and hence by the Stone-Weierstrass theorem,

[TABLE]

Since $F$ and $C$ are nonvanishing, equation (51) implies that $\dot{\rho}=0$.

6.2. Proof of $\dot{K}=0$

We now return to the first variation of the Robin function $K$. According to Theorem 3.7,

[TABLE]

Since we have just shown that $\dot{\rho}=0$ in Section 6.1, the variational trace formula in Theorem 3.7 becomes

[TABLE]

In Section 5.3, we cooked up an oscillatory integral representation for $S_{j}$ microlocally near the canonical relations $\Gamma_{\pm}^{j}$ lying over an open neighborhood of the diagonal of the boundary (Theorem 5.7). Without the $L_{R}^{b}\dot{\rho}$ terms, we do not differentiate the sine kernel in the integrand of (52). Again following the formulas for homogeneous distributions in [Hör03], we see that near $t=T_{j}$,

[TABLE]

where the minus sign is due to the appearance of the same Maslov index $\sigma$ appearing in the principal term. Plugging in the formula (44) for $A_{j}$ gives the second formula in Theorem 1.4. Similar computations to those in Section 6.1 lead us to the equation

[TABLE]

Discarding the separated terms depending only on $\lambda_{j}$ and Taylor expandng at $\lambda=0$ as before, we see via the Stone-Weierstrass theorem that $\dot{K}=0$, which concludes the proof of Theorem 1.1.

Remark 6.1*.*

In [HZ12], the principal symbol computation for Neumann boundary conditions follows closely the work of [DG75]. The principal symbols for the Neumann and Robin wave propagators agree and the computations in [HZ12] yield

[TABLE]

where the sum is over connected components $\Gamma$ of the fixed point set $F_{T}$ and $\gamma_{1}\in C(B^{*}\partial\Omega)$ is given by $\gamma_{1}(q,\zeta)=\sqrt{1-|\zeta|^{2}}$. The coefficients $C_{\Gamma}$ are nonzero Maslov factors coming from the stationary phase computation in [HZ12] and $\mu_{L}$ is the Leray measure on $\Gamma\subset B^{*}\partial\Omega$, which is computed in [GM79a]. In our computations, all of these factors are explicit.

7. Proof of Lemma 5.3

This section is dedicated to the proof of Lemma 5.3. Fix $j\in\mathbb{Z}$ large and choose a corresponding tubular neighborhood $U$ of $\partial\Omega$ with the following property: for any $x\in U$, the $j$-fold broken geodesic emanating from $x$ which is tangent to the confocal ellipse on which $x$ lies makes less than a quarter rotation. By this we mean that if $x=(a\cosh\mu_{x}\cos\varphi_{x},b\sinh\mu_{x}\sin\varphi_{x})$, then the $j+1$ impact point at the boundary has angular component in the interval $\varphi_{x}\leq\varphi_{x,j}\leq\varphi_{x}+\pi/2$. This is certainly possible since the billiard flow is continuous on the closure of its phase space and if $x$ is on the boundary, the corresponding orbit is stationary. Denote by $C_{\lambda_{x}}$ the confocal ellipse on which $x$ lies. We will perturb this orbit by holding $x$ fixed and increasing the parameter $\lambda^{2}$ of the confocal caustic

[TABLE]

to which the orbit is tangent. This can be done by rotating the initial covector of the trajectory slightly in any direction within $S_{x}^{*}(\mathbb{R}^{2})$ so that it makes a nonzero angle with the tangent line to $C_{\lambda_{x}}$ at $x$. As $\lambda^{2}$ increases, the associated angle also increases, the confocal ellipses shrink and heuristically, the $j$-fold broken geodesic begins to rotate more and more around $\Omega$. This is precisely the twist property of the billiard map: for a fixed point $x$ in the base, the straight line in phase space obtained by letting the angular component vary becomes twisted under iteration of the billiard map. If $y$ is another point which is sufficiently close to $x$, we claim that there exist $8$ angles in $S_{x}^{*}(\mathbb{R}^{2})$ such that the last link of the corresponding billiard trajectory intersects $y$.

Four of these orbits will be oriented in the counterclockwise direction and four in the clockwise direction. Of the four counterclockwise orbits, two of them will correspond to rotating the initial covector in $S_{x}^{*}(\mathbb{R}^{2})$ in the counterclockwise direction and two will result from a rotation in the opposite direction. We show that in elliptical polar coordinates, the angular components of both intersection points of the last link (after $j$ reflections) with the confocal ellipse on which $y$ lies are increasing as we rotate the inital covector within either direction in $S_{x}^{*}(\mathbb{R}^{2})$. Hence, both intersection points will wind around the confocal ellipse until they eventually coincide with $y$. The clockwise orbits will then be constructed by a simple reflection argument.

7.1. Notation

Let $\lambda_{x}$ and $\lambda_{y}$ denote the parameters of the confocal ellipses on which $x$ and $y$ lie respectively. If $\lambda_{x}>\lambda_{y}$, then it is clear that the last link of any billiard emanating from $x$ will intersect $C_{\lambda_{y}}$ exactly twice. If $\lambda_{y}>\lambda_{x}$ but $x$ and $y$ are sufficiently close, then $\lambda_{x}$ and $\lambda_{y}$ are also close, so we can arrange that the billiard emanating from $x$ which is tangent to $C_{\lambda_{y}}$ makes less than a half rotation. Hence, the orbits making approximately one full rotation will necessarily have final links which intersect $C_{\lambda_{y}}$ twice. We only consider such $x,y$, which lie in an open neighborhood of $\Delta\partial\Omega\subset\Omega\times\Omega$, where $\Delta:\partial\Omega\to\partial\Omega\times\partial\Omega$ is the diagonal embedding.

The aim is to prove that the angular components of both intersection points of the last link with the caustic $C_{\lambda_{y}}$ are increasing. To do this, we will consider a variant of the map $B_{j}$ defined in Section 5.3. Recall that we defined

[TABLE]

where $\beta$ is the billiard map on the coball bundle $B^{*}(\partial\Omega)$ of the boundary. Since the coball bundle of the boundary is diffeomorphic to the collection of inward facing covectors in $S_{\partial\Omega}^{*}(\mathbb{R}^{2})$ and the latter is more geometrically natural, we lose no generality by considering $B_{j}$ or $\beta$ as a map on $S_{\partial\Omega}^{*}(\mathbb{R}^{2})$. We’ve assumed that $x,y\in\text{int}\Omega$, so we must first flow to the boundary in order to study the billiard map. Let $\alpha\in S^{1}$ and denote by $(x_{1}(\alpha),\omega_{1}(\alpha))$ the point obtained by evolving $(x,\alpha)$ under the forward billiard flow for $t_{+}^{1}(x,\alpha)$ units of time and then reflecting at the boundary. Recall from Section 4 that we defined

[TABLE]

To the point $(x_{1}(\alpha),\omega_{1}(\alpha))\in S_{\partial\Omega}^{*}(\mathbb{R}^{2})$ which is now fixed, we may apply $B_{j-1}$ to obtain the angular component of the $j$th impact point. Denote the corresponding boundary point by

[TABLE]

Denote by $\alpha_{x}$ the angle which the positively oriented tangent line $T_{x}C_{\lambda_{x}}$ makes with the positive $x$ axis. If $\alpha$ is sufficiently close to $\alpha_{x}$, the first $j$ iterates under billiard map will also make less than a quarter rotation by our original set up. Our strategy is to let $\alpha$ vary within an open cone of directions which parametrizes all possible counterclockwise orbits emanating from $x$, and show that precisely $4$ of these angles result in orbits which reach $y$ after $j$ reflections, making approximately one rotation.

For each $\alpha$, there exists a unique $\lambda$ such that the corresponding orbit is tangent to the caustic $C_{\lambda}$. For example, $\lambda_{x}$ corresponds to $\alpha_{x}$.

Definition 7.1**.**

Let $U_{x}^{*}$ be the set of $\alpha\in S_{x}^{*}(\mathbb{R}^{2})$ such that the corresponding orbits avoid the region between the focal points. Then, denote the homogeneous extension of $U_{x}^{*}$ to $T_{x}^{*}\mathbb{R}^{2}$ by $C_{x}^{*}$, which is precisely the fiber cone at $x$ of admissible initial covectors we consider.

Remark 7.2*.*

If $x$ varies smoothly in $\Omega$, so do the associated fibers and in this way, we obtain a smooth cone bundle, which we denote by $C^{*}(\Omega)\subset T^{*}(\Omega)$, over a tubular neighborhood of $\partial\Omega$.

7.2. Derivatives of $\beta$

If $\alpha$ increases, so does the parameter $\lambda$ of the confocal caustic, i.e. $d\lambda/d\alpha>0$ for $\alpha>\alpha_{x}$. Similarly, as $\alpha$ decreases, the parameter $\lambda$ also increases, as can be seen by considering the backwards (clockwise) orbit. Hence, $\lambda_{x}$ is a local minimum for $\alpha$ near $\alpha_{x}$. By the implicit function theorem, it is clear that $\lambda$ is a smooth function of $\alpha$ as long as the corresponding forwards and backwards orbits do not enter between the focal points. Let us first consider the case in which $\alpha$ increases. For $\alpha>\alpha_{x}$, there is a one to one correspondence between $\alpha$ and $\lambda$. In this case, it is geometrically clear that the angular component of $x_{1}$ is also increasing in $\alpha$ and $\partial x_{1}/\partial\alpha$ is bounded independently of $j,\lambda$, $\omega$.

To obtain a large lower bound on the speed at which the two intersection points wind around $C_{\lambda_{y}}$, we first find a corresponding lower bound on

[TABLE]

By a slight abuse of notation, we have systematically confused $x_{1}$ with its angular component $\theta(x_{1})$, since they are in one to one correspondence modulo factors of $2\pi$. We begin with a simple lemma which will be used throughout the section.

Lemma 7.3**.**

Let $x,y\in\Omega$ be $O(1/j)$ close to the diagonal of the boundary and $\gamma$ be an orbit of $j$ reflections which is tangent to the caustic $C_{\lambda}$ and connects $x$ to $y$. Also denote by $\omega_{k}$ ($1\leq k\leq j$) the angle of reflection made at the $k$th impact point on the boundary. Then $\lambda=O(1/j)$ and $\omega_{k}=O(1/j)$ for all $1\leq k\leq j$.

Proof.

Recall that by $\gamma$ making approximately one rotation we mean that $|\widehat{\beta}(\widehat{x},\omega_{1})-\widehat{y}-\ell|<\ell/100$ (see the remarks following Lemma 5.3). If in action angle coordinates (cf. Section 5.5), $\widehat{x}$ is given by $q_{\lambda}(s)=(-a\text{sn}(s,k_{\lambda}),b\text{cn}(s,k_{\lambda}))$ for some $s\in\mathbb{R}$, we have that the link

[TABLE]

is tangent to the confocal ellipse $C_{\lambda}$. Lifting to the universal cover $\mathbb{R}\times(0,\pi)$ as in Section 4 and setting $\widetilde{s}\in\mathbb{R}/\ell\mathbb{Z}$ to be the arclength parameter corresponding to $s$, we have $\pi_{1}\widehat{\beta}(\widetilde{s},\omega_{1})=\widetilde{s}+\delta_{\lambda}$ where $\pi_{1}$ is projection onto the first factor. Similarly,

[TABLE]

If $y$ is given by $(-a\text{sn}(t,k_{\lambda}),b\text{cn}(t,k_{\lambda}))$ in action angle coordinates and $\widetilde{t}$ is the associated arclength parameter, then $x$ being $O(1/j)$ close to $y$ and the boundary is equivalent to requiring $|\widetilde{s}-\widetilde{t}|=O(1/j)$. Then (54) implies that

[TABLE]

Hence, $\delta_{\lambda}=\ell/j+O(j^{-2})$. Note that

[TABLE]

As ${\sqrt{1-\lambda^{2}k_{\lambda}^{2}/b}}=1+O(\lambda)$ and $\arcsin\lambda/b=\lambda/b+O(\lambda^{2})$, $\delta_{\lambda}=O(1/j)$ immediately implies that $\lambda=O(1/j)$. The relationship between $\omega_{1}$ and $\lambda$ is given by Lemma 5.14: $\lambda^{2}=\sin^{2}\omega_{1}(b^{2}+(a^{2}-b^{2})\sin^{2}\widehat{x})$. As the coefficient of $\sin^{2}\omega_{1}$ is bounded above by $a$ and below by $b$, $\lambda$ and $\omega_{1}$ are of the same order near zero. Furthermore, each link of $\gamma$ is tangent to $C_{\lambda}$ for a fixed $\lambda=O(1/j)$ and hence, the same logic also implies that $\omega_{k}=O(1/j)$ for each $1\leq k\leq j$. This concludes the proof of the lemma. ∎

We are now ready to estimate (53).

Lemma 7.4**.**

For $j$ large, we have

[TABLE]

Proof.

In Section 5.3, we found a rather explicit expression for

[TABLE]

where $x_{1}\in\partial\Omega$ is fixed and $\omega$ is evaluated at the critical angle ${\omega}(x_{1})$ corresponding to a periodic orbit. Using the formulas in Section 5.3, we can actually let $x_{1}$ vary in $\alpha$ and evaluate at any $\omega$; in particular, we can evaluate at $\omega=\omega_{1}$ corresponding to the angle of reflection at the first impact point on the boundary. Let us first examine the term in (57) with a coefficient of $(j+1)$. We see that

[TABLE]

If $j$ is suffciently large, we can ensure that $\lambda=O(1/j)$ is in turn small. Hence, (58) can bounded below by $(1-b^{2}/a^{2})/2$ and recalling Lemma 5.14, the product of $1/\lambda$ and $\sin\omega$ coming from equations (59) and (60) can be estimated below by $1/a$. For $\lambda\ll 1$, all of the remaining terms can easily be bounded below by a positive constant depending only on $a$ and $b$. Hence, we have

[TABLE]

We now consider the first term $\frac{\partial\text{am}}{\partial s}\frac{\partial t_{\varphi}}{\partial\lambda^{2}}\frac{\partial\lambda^{2}}{\partial\omega}$ in (57). The first and third factors of this product are clearly bounded above and below, independently of $\omega,\lambda$, and $j$ near $\lambda=0$ by the same arguments as before. Recall from Section 5.3 that

[TABLE]

which is also clearly bounded in magnitude by a positive constant independent of $\omega,\lambda$ and $j$.

In a similar manner, we would like to estimate the final term $\frac{\partial\text{am}}{\partial k_{\lambda}^{2}}\frac{dk_{\lambda}^{2}}{d\lambda^{2}}\frac{\partial\lambda^{2}}{\partial\omega}$ in (57). Recall from Section 5.3 that

[TABLE]

Both terms can be bounded independently of $\omega,\lambda$ and $j$. Combining this with the earlier bound for (60), we see that

[TABLE]

We also need to estimate the remaining terms in (53). In particular, we must bound

[TABLE]

from below by a positive constant and

[TABLE]

in magnitude. Lemma 5.14 tells us that if $0\leq\omega\leq\pi/3$, then

[TABLE]

so that

[TABLE]

It is also geometrically clear that

[TABLE]

By writing out $B_{j}$ in action angle coordinates, we see that

[TABLE]

where $t_{\varphi}$ is defined implicitly by the equation $\varphi=\text{am}(t_{\varphi};k_{\lambda})+\pi/2$. Hence, by the implicit function theorem and our previous bound for (58), $\frac{\partial B_{j}}{\partial x_{1}}$ is both nonnegative and bounded independently of $\omega,\lambda$ and $j$. Combining (61), (62), (63) and (64), we see that

[TABLE]

∎

7.3. Intersection points

We are now ready to prove that the intersection points with the last link wind monotonically around the caustic $C_{\lambda_{y}}$ on which $y$ lies as we increase $\lambda$. In tandem with Lemma 7.4, we will also need information about how the final angle $\omega_{j}$ at the $j$th reflection depends on $\lambda$.

Lemma 7.5**.**

If $\omega_{j}=\pi_{2}\beta^{j}(x_{1},\omega_{1})$ denotes the angle of reflection at the $j$th impact point on the boundary, then $\partial\omega_{j}/\partial\lambda=O(1)$, i.e. it’s derivative bounded independently of $j$.

Proof.

Recall Lemma 5.14, which gave

[TABLE]

Differentiating this as in the proof of Lemma 7.4 but now using that $x_{j}$ depends on $\lambda$, we see that

[TABLE]

Substituting again the formula in Lemma 5.14 into equation (65), we find that

[TABLE]

Lemma 7.3 tells us that $\lambda,\omega_{j}=O(1/j)$ and hence $\cos\omega_{j}=1+O(1/j),\sin\omega_{j}=O(1/j)$. This implies that the denominator of (66) is bounded below by a positive constant for $j$ larger than some fixed $j_{0}$ depending only on $\Omega$. The only unbounded term in the numerator is $\partial x_{j}/\partial\lambda=\partial B_{j}/\partial\lambda\gtrsim j$, as was shown in Lemma 7.4. However, an examination of the proof of Lemma 7.4 actually shows that $\partial x_{j}/\partial\lambda=O(1/j)$ and hence, $\sin\omega_{j}\frac{\partial x_{j}}{\partial\lambda}=O(1)$, which implies that (66) is in fact bounded. ∎

With Lemmas 7.3, 7.4 and 7.5, we can prove half of Lemma 5.3:

Lemma 7.6**.**

For $j$ sufficiently large, the angular components in elliptical polar coordinates of both intersection points $\varphi_{j}^{1}$ and $\varphi_{j}^{2}$ of the final link with the caustic $C_{\lambda_{y}}$ are monotonically increasing in $\lambda$ with approximate speed $j$.

Proof.

Recall from the proof of Lemma 5.14 that the billiard ray emanating from $x_{j}$ at angle $\omega_{j}$ is parametrized by

[TABLE]

Taking real and imaginary parts, we find that

[TABLE]

For simplicity, denote by $A$ the coefficient of $t$ in (67). Converting the line $L$ from parametric form, we find that

[TABLE]

We know that $L$ has exactly two intersection points with

[TABLE]

the caustic on which $y$ lies. These correspond to two different values of $t$ in equation (67). At either intersection point, the values of $x,y$ in (68) are constrained to be of the form $(c\cosh\mu_{j}\cos\varphi_{j},c\sinh\mu_{j}\sin\varphi_{j})$, where $c=\sqrt{a^{2}-b^{2}}$, $\mu_{j}=\cosh^{-1}((a^{2}-\lambda_{y}^{2})^{1/2}/c)$ and $\varphi_{j}\in[0,2\pi)$ is the angular parameter in elliptical polar coordinates for the caustic $C_{\lambda_{y}}$. For simplicity, denote

[TABLE]

To begin, assume that $\varphi_{j}$ is either the first intersection point or the second intersection point. Only at the end of the proof will we need to distinguish between the two cases. We want to show precisely that $\varphi_{j}^{\prime}=\frac{\partial}{\partial\lambda}\varphi_{j}>0$ so that the intersection points wind around $C_{\lambda_{y}}$. Solving for $w$ in equation (68) and substituting $z=c_{j}\cos\varphi_{j},w=s_{j}\sin\varphi_{j}$, we see that

[TABLE]

Differentiating (69) and collecting terms, we have

[TABLE]

where

[TABLE]

and

[TABLE]

As both intersection points are within $O(1/j)$ of $x_{j}$, we have also that $|x_{j}-\varphi_{j}|=O(1/j)$ and $D_{3}=O(1/j)$. However, setting $\omega_{j}=0$, $\varphi_{j}=x_{j}$, $c_{j}=a$, and $s_{j}=b$, we see that $D_{1}$, $D_{2}$ and $D_{3}$ all vanish. Instead we Taylor expand each coefficient to second order in $1/j$. Recall that we are free to choose $y$ as close to the boundary as we want, so $\lambda_{y}=O(j^{-N})$, $c_{j}=a+O(j^{-N}),$ and $s_{j}=b+O(j^{-N})$ for any $N\in\mathbb{N}$. Expanding $D_{1}$ and simplifying, we see that

[TABLE]

Similarly, we obtain

[TABLE]

and

[TABLE]

Hence, $D_{1}=D_{2}+O(j^{-2})$ and $D_{3}=O(1/j)$. As $x_{j}^{\prime}\gtrsim j$ by Lemma 7.4, we are done as long as we can show that $D_{1},D_{2}$ don’t vanish to first order in $1/j$ so that we may divide through by them in equation (70). If $\varphi_{j}$ is the first intersection point, then $\varphi_{j}=x_{j}+O(j^{-N})$ for any $N$ so modulo $O(j^{-2})$ terms, $D_{1}$ becomes $-\omega_{j}(a^{2}\sin^{2}x_{j}\varphi+b^{2}\cos^{2}x_{j})\leq-c\omega_{j}\leq-c^{\prime}/j$ for some constants $c,c^{\prime}>0$ depending only on the curvature of $\partial\Omega$. Hence, we may assume $\varphi_{j}$ is the second intersection point, in which case choosing $y$ sufficiently close to $\partial\Omega$ amounts to setting $\varphi_{j}=x_{j+1}+O(j^{-N})$ for any desired $N$.

Simplifying further, we see that modulo terms of order $O(j^{-2})$,

[TABLE]

We want to show that the positivity of the first term in (74) outweights the second term. Recall from Section 5.5 that $x_{j+1}-x_{j}=\text{am}(s+\delta_{\lambda},k_{\lambda})-\text{am}(s,k_{\lambda})$ for some $s\in\mathbb{R}$ and $\lambda$ (not $\lambda_{y}$) corresponding to the caustic to which the orbit is tangent. Hence,

[TABLE]

which implies that

[TABLE]

As $\lambda=O(1/j)$ for $j$ large, $\arcsin\lambda/b=\lambda/b+O(j^{-2})$ and $(1-\lambda^{2}k_{\lambda}^{2}/b)^{1/2}=1+O(\lambda)=1+O(1/j)$. Hence, modulo terms of order $O(j^{-2})$, we have

[TABLE]

Now recall Lemma 5.14, which gives

[TABLE]

Near $\omega_{j}=0,\lambda=0$, this tells us that

[TABLE]

The term in $D_{1}$ with a $\sin\omega_{j}$ can hence be estimated modulo $O(j^{-2})$ by

[TABLE]

The coefficient of $\lambda$ in (78) is positive and maximized at critical points corresponding to $\sin x_{j}=0,\cos x_{j}=0$ or

[TABLE]

If $\sin x_{j}=0$, then $\alpha_{j}=\lambda/b+O(j^{-2})$ and the leading order coefficient of $\omega_{j}$ in the formula for $D_{1}$ is $b^{2}$. Hence, (77) implies that

[TABLE]

and $D_{1}>0$ is nonvanishing. If $\cos x_{j}=0$, then $\omega_{j}=\lambda/a+O(j^{-2})$ and the leading order coefficient of $\omega_{j}$ in the formula for $D_{1}$ is $a^{2}$. Hence, (77) again implies that

[TABLE]

and $D_{1}>0$ is nonvanishing. In the third case, solving (79) results in

[TABLE]

which impossible unless $b^{2}=0$. ∎

In this way, we obtain two orbits from $x$ to $y$ by letting $\alpha$ increase within the $\alpha\geq\alpha_{x}$ regime. Dynamically, these orbits can be characterized by having a point of tangency to a confocal ellipse before a moment of reflection at the boundary in the forwards direction.

We now consider the regime in which $\alpha\leq\alpha_{x}$. There is a different one to one correspondence between $\alpha$ and $\lambda$ but many of the equations above remain completely valid. In particular (61) and (62) are unchanged. Hence, by essentially the same bounds as in Lemma 7.4 before, if $\varphi_{j}^{(1)}$ and $\varphi_{j}^{(2)}$ denote the angular components of the intersection points in elliptical polar coordinates, we have

[TABLE]

in the $\alpha\leq\alpha_{x}$ regime. This provides two additional orbits, which are dynamically characterized by having a reflection at the boundary before becoming tangent to a confocal ellipse.

To obtain the four clockwise orbits, note that we can first apply the isometry $\mathbb{R}^{2}\ni(z,w)\mapsto(-z,w)$, obtain four counterclockwise orbits as above, and then reflect back. In the last link of any orbit, the first point of intersection with $C_{\lambda_{y}}$ is reached before a point of tangency with a confocal ellipse while the second intersection point is reached after a point of tangency. These characterizations are important in Section 5.3 for understanding four different types of orbit configurations and determining which types of limiting orbits give rise to periodic orbits of precisely $j$ reflections as $(x,y)$ approach the diagonal of the boundary.

Remark 7.7*.*

For sufficiently large $j$ and $x,y$ both lying on the boundary near the diagonal, the existence of a single such geodesic for general smooth, strictly convex domains was proven in [GM81], [MM82] and [Pop94]. The eight orbits in the statement of Lemma 5.3 can be seen to collapse into the orbits described in [MM82] as $(x,y)$ approach the diagonal of the boundary from any direction. However, there may be a different number of reflections in the limiting orbit (see the proof of Lemma 5.8).

8. Acknowledgements

The author would like to thank Hamid Hezari and Katya Krupchyk for their support and suggestions during this project. The author would also like to thank Vadim Kaloshin and Alfonso Sorrentino for allowing the use of their images in Figure 1. Additionally, Katya Krupchyk was gracious enough to provide the author with funding from NSF grant DMS 1500703 during the summer of 2016.