Balian-Bloch wave invariants for nearly degenerate orbits

Vadim Kaloshin ¹ , Illya Koval ¹ and Amir Vig ² ¹ Institute of Science and Technology Austria, Klosterneuburg, Lower Austria ² Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA Vadim.Kaloshin@gmail.com Illya.Koval@ist.ac.at Vig@umich.edu

Abstract.

This paper is part I of a series in which we aim to show that the singular support of the wave trace and the length spectrum of a smooth, strictly convex, and bounded planar billiard table are generally distinct objects. We derive an asymptotic trace formula for the regularized resolvent which is dual to the wave trace and contains the same information. To do this, we consider a class of periodic orbits which have nearly degenerate Poincaré maps and study their leading order behavior as the deformation parameter goes to zero, generating large coefficients in the wave trace. We also keep careful track of the Maslov indices, which will allow us to match contributions of opposite signs in our subsequent paper. Each cancellation of coefficients in the resolvent trace corresponds to making the wave trace one degree smoother. The resolvent based approach is due to Balian and Bloch and was significantly expanded upon by Zelditch in a foundational series of papers [Zel09], [Zel04a], [Zel04c] and [Zel00].

1. Introduction

In this paper, we consider the Laplacian on a smooth, strictly convex and bounded planar domain $\Omega$ :

(1)

\displaystyle\begin{cases}-\Delta u=\lambda^{2}u&x\in\Omega\\ Bu=0,&x\in\partial\Omega,\end{cases}

$B$ is a boundary operator corresponding to Dirichlet, Neumann or Robin boundary condtions. The even wave trace $w(t)$ is defined by the real part of the distribution

(2)

\displaystyle\frac{1}{2\pi}\int_{-\infty}^{+\infty}\text{Tr}e^{it\sqrt{-\Delta}}\hat{\varphi}(t)dt=\text{Tr}\varphi(\sqrt{-\Delta}),\qquad\varphi\in\mathcal{S}(\mathbb{R})

and we write $w(t)=\text{Tr}\cos({t\sqrt{-\Delta}})$ . Denote by $\text{LSP}(\Omega)$ the (unmarked) length spectrum, consisting of lengths of periodic billiard trajectories in $\Omega$ . The Poisson relation tells us that the singular support of $w(t)$ is contained in the closure of the length spectrum together with $0$ and its reflection about the origin. Assume $L$ is isolated in $\text{LSP}(\Omega)$ with finite multiplicity and all corresponding periodic orbits are nondegenerate. Choose $\widehat{\rho}\in C_{c}^{\infty}$ to be identically equal to $1$ in a neighborhood of $L$ and satisfy $\text{supp}(\widehat{\rho})\cap\text{LSP}(\Omega)=L$ . We then have an asymptotic expansion of the regularized resolvent trace which has form

(3)

\displaystyle\int_{0}^{\infty}\widehat{\rho}(t)w(t)e^{ikt}dt\sim\sum_{\gamma:\text{length}(\gamma)=L}\mathcal{D}_{\gamma}\sum_{j=0}^{\infty}B_{\gamma,j}k^{-j},

where the righthand side is a sum over all periodic orbits $\gamma$ having length $L$ . We call $\mathcal{D}_{\gamma}$ the symplectic prefactor and $B_{\gamma,j}$ the Balian-Bloch wave invariants.

Theorem 1.1.

Let $\Omega_{\varepsilon}$ , $\varepsilon\in[0,\varepsilon_{0}]$ be a smooth one parameter family of domains which fixes the reflection points and angles of a period $q\geq 3$ billiard orbit $\gamma$ . Assume further that $\partial^{2}\mathcal{L}$ has rank $q-1$ in $\Omega_{0}$ , where $\partial^{2}\mathcal{L}$ is the Hessian of the length functional in arclength coordinates evaluated at $\gamma$ , and that the perturbation makes $\gamma$ nondegenerate in $\Omega_{\varepsilon}$ when $\varepsilon>0$ , with $|\det\partial^{2}\mathcal{L}|\sim c_{\gamma}\varepsilon$ for some $c_{\gamma}>0$ . Then, the $j$ th Balian-Bloch wave invariant $B_{j,\gamma}$ associated to $\gamma$ in the regularized resolvent trace 3 is given by

\displaystyle\frac{e^{-\frac{i\pi q}{4}}}{2q}L_{\gamma}\left(\sum_{i_{1},i_{2},i_{3}=1}^{q}h_{i_{1}}h_{i_{2}}h_{i_{3}}\partial_{i_{1},i_{2},i_{3}}^{3}\mathcal{L}\right)^{2j}(c_{\gamma}\varepsilon)^{-3j}\prod_{\ell=1}^{q}\frac{\cos\vartheta_{\ell}}{|x(s_{\ell})-x(s_{\ell+1})|}\frac{(\pm i)^{j}}{w(j)}+\mathcal{R}_{j},

where

•

$|x(s_{\ell})-x(s_{\ell+1})|$ is the length of the $\ell$ th link in the billiard trajectory.
•

$\vartheta_{\ell}$ are the angles of reflection measured from the interior normal at $\ell$ th impact point on the boundary.
•

$w(j)=\sum_{\mathcal{G}}\frac{1}{w(\mathcal{G})}>0$ is a sum is over all $3$ -regular graphs on $2j$ vertices, with $w(\mathcal{G})$ being the order of the automorphism group of a graph $\mathcal{G}$ . In particular, it is independent of both the domain and the orbit.
•

$\mathcal{L}$ is the length functional in arclength coordinates (see Definition 3.1), $L_{\gamma}$ is the length of $\gamma$ and the functions $h_{i}\in C^{\infty}(\partial\Omega_{0})$ are such that $h_{i_{1}}h_{i_{2}}$ is, up to a sign, the $(i_{1},i_{2})$ entry of the adjugate matrix of $\partial^{2}\mathcal{L}$ at $\gamma$ in $\Omega_{0}$ , which is independent of the deformation. The factor

$\sum_{i_{1},i_{2},i_{3}}h_{i_{1}}h_{i_{2}}h_{i_{3}}\partial_{i_{1},i_{2},i_{3}}^{3}\mathcal{L}$

is, modulo an $O(\varepsilon)$ error, the third directional derivative of $\mathcal{L}$ along the null space of $\partial^{2}\mathcal{L}$ in $\Omega_{0}$ .
•

The signs $\pm$ are given by $+$ if the closest eigenvalue to $0$ of $\partial^{2}\mathcal{L}$ in $\Omega_{\varepsilon}$ is positive and $-$ otherwise.
•
$\mathcal{R}_{j}=\mathcal{R}(\gamma,\varepsilon,j,\kappa,\kappa^{(1)},\cdots,\kappa^{(2j)})$ , with $\kappa^{(i)}$ being the $i$ th derivative of the boundary curvature of $\Omega_{\varepsilon}$ , is a remainder satisfying:
- –
  
  $\mathcal{R}_{j}=O(\varepsilon^{-3j+1})$ .
- –
  
  $\varepsilon^{3j-1}\mathcal{R}_{j}$ is smooth in all parameters down to $\varepsilon=0$ .
- –
  
  $\mathcal{R}_{j}$ depends on at most $2j$ derivatives of the boundary curvature in a small neighborhood of the reflection points of $\gamma$ .
•

When $j=0$ , $w(0)=1$ and $\mathcal{R}_{0}=0$ .

The symplectic prefactor associated to $\gamma$ is

\displaystyle\mathcal{D}_{\gamma}(k)=\frac{e^{ikL_{\gamma}}e^{\frac{i\pi\text{sgn \,}\partial^{2}\mathcal{L}(S_{\gamma})}{4}}}{\sqrt{|\det\partial^{2}\mathcal{L}(S_{\gamma})|}}.

Our paper is part of a series in which we aim to show that the singular support of the wave trace and the length spectrum are in general distinct objects which encode different information. In part II, we will create cancellations in the wave trace by perturbing an initial domain having the properties in Theorem 1.1 and then using the asymptotic formulas above together with a careful matching of Maslov indices, which arise as complex phases. The possibility of a smooth wave trace was first posed in [DG75] and it was later remarked by Colin de Verdiere (see [Zel04b]) that apriori, there is no immediate reason why $w(t)$ cannot be $C^{\infty}$ on $\mathbb{R}\backslash\{0\}$ ; it is known to be singular at $t=0$ . In a work in progress, Hezari and Zelditch construct domains for which pairs of bouncing ball orbits give infinite order cancellations in the wave trace, going further in this direction. Our result applies to a large class of convex domains, having a rather general class of degenerate orbits; see Section 5.

2. Background

Historically, the wave trace invariants have been exploited in the context of the inverse spectral problem, first posed by Mark Kac in [Kac66]: “Can one hear the shape of a drum?” Mathematically, this asks if one can recover the geometry of a domain from knowledge of its Laplace spectrum. In general the answer is no, as was shown in [Mil64], [GWW92], [Vig80], and [SM85] amongst others. However, there are many cases in which one expects the answer to be yes, perhaps the most tantalizing of which is the class of strictly convex, smooth planar domains, which we consider in the present article. There have been numerous works in this direction, perhaps the most significant of which are [Zel09] and [HZ22]. In the first article, Zelditch showed that generic families of analytic $\mathbb{Z}_{2}$ symmetric planar domains can be recovered from their spectrum by a careful analysis of the wave invariants associated to iterates of a bouncing ball orbit. We follow parts of this approach closely in our paper. In the second, Hezari and Zelditch prove that ellipses of small eccentricity are spectrally determined amongst all smooth, strictly convex bounded domains by their spectrum. This is the first and only known example of noncircular domains which are spectrally determined. $C^{\infty}$ compactness of Laplace isospectral sets was shown in [OPS88c], [OPS88b] and [OPS88a]. Other recent progress can be found in [Zel00], [Zel98], [Zel04c], [HZ10], [HZ12], [Vig21], and [Vig23b].

The dynamical analogue of the inverse spectral problem is to determine the shape of a manifold with boundary from knowledge of the length spectrum, consisting of lengths of closed geodesics (billiard trajectories) which make elastic reflections at the boundary. For planar domains, recent progress in this direction was acheived by the first author and Sorrentino in [KS18], together with the works [ADSK16], [Kov21] and [dSKW17]. In the first three, it was shown that ellipses are isolated within the class of integrable domains. The latter two concern spectral rigidity of ellipses and nearly circular domains, respectively. $C^{\infty}$ compactness of marked length isospectral sets was demonstrated by the third author in [Vig23a]. Other significant results, including cases without boundary, are contained in [HKS18], [PT12], [Pop94], [PT11], [GK80], [CdV84] and [GL18].

3. Billiards

Before obtaining a singularity expansion for the wave trace, we first review the relevant background on billiards. This will be useful in our construction of the Balian-Bloch resolvent parametrix in Section 4.2. We denote by $\Omega$ a bounded and strictly convex region in $\mathbb{R}^{2}$ with smooth boundary. This means that the curvature of $\partial\Omega$ is a strictly positive function. We can identify $B^{*}\partial\Omega$ with ${\mathbb{R}}/{\ell\mathbb{Z}}\times(0,\pi)$ , where $\ell=|\partial\Omega|$ is the length of the boundary. Let $x:[0,\ell]\to\mathbb{R}^{2}$ be a unit speed parametrization of $\partial\Omega$ and denote by $ds$ the arclength measure on $\partial\Omega$ . The billiard map is defined on the coball bundle of the boundary, $B^{*}\partial\Omega=\{s,\sigma)\in T^{*}\partial\Omega:|\zeta|<1\}$ , which can be identified in two ways with the cosphere bundle $S_{\partial\Omega}^{*}\mathbb{R}{{}^{2}}$ over the boundary via the two orthogonal projection maps from inward $(+)$ /outward $(-)$ pointing covectors. For $\sigma\in B^{*}\partial\Omega$ , define $\varphi_{\pm}$ to be the inward $(+)$ /outward $(-)$ facing covectors in $S_{\partial\Omega}^{*}\mathbb{R}^{2}$ which project to $\sigma$ .

Define the maps

\displaystyle\beta_{\pm}(s,\sigma)=(s^{\prime},\sigma^{\prime}),

where $s^{\prime}$ is the arclength coordinate of the subsequent intersection of the line through $x(s)$ with direction $\varphi_{\pm}$ and $\sigma$ is the projection onto $B^{*}\partial\Omega$ of the parallel transport of $\varphi_{\pm}$ to $x(s^{\prime})$ . We call $\beta=\beta_{+}$ the billiard map. It is well known that $\beta$ preserves the canonical one form $\sigma ds$ induced on $B^{*}(\partial\Omega)$ and is differentiable there, extending continuously up to the boundary, with a square root type singularity at $0$ and $\pi$ corresponding to a Whitney fold in its graph. As a consequence, $\beta$ also preserves the area form $d\sigma\wedge ds$ . The maps $\beta_{\pm}^{n}$ are defined via iteration and it is clear that $\beta_{\pm}^{-n}=\beta_{\mp}^{n}$ for $n\in\mathbb{Z}$ . It is also apparent that $\frac{\partial\beta}{\partial\vartheta}>0$ , which means that vertical fibers in $B^{*}\partial\Omega$ are twisted upon iterating $\beta$ . Together, these properties make $\beta$ into an exact symplectic twist map. It has a particularly simple generating function

\displaystyle h(s,s^{\prime})=-|x(s)-x(s^{\prime})|,

by which mean that $\beta\left(x,-\frac{\partial h}{\partial s}\right)=\left(s,\frac{\partial h}{\partial s^{\prime}}\right)$ , which implies that $s$ and $s^{\prime}$ give local coordinates on the graph of $\beta$ . Geometrically, a billiard orbit corresponds to a union of line segments which are called links. A smooth closed curve $\mathcal{C}$ lying in $\Omega$ is called a caustic if any link drawn tangent to $\mathcal{C}$ remains tangent to $\mathcal{C}$ after an elastic reflection at the boundary of $\Omega$ . By elastic reflection, we mean that the angle of incidence equals the angle of reflection at an impact point on the boundary. We can map $\mathcal{C}$ onto phase space $B^{*}\partial\Omega$ to obtain a smooth closed curve which is invariant under $\beta_{\pm}$ . If the dynamics are integrable, these invariant curves are precisely the Lagrangian tori which folliate phase space. A point $P$ in $B^{*}\partial\Omega$ is called $q$ -periodic ( $q\geq 2$ ) if $\beta^{q}(P)=P$ . We define the rotation number of a $q$ -periodic orbit $\gamma$ emanating from $P\in\overline{B^{*}(\partial\Omega)}$ by $\omega(\gamma)=\omega(P)=\frac{p}{q}$ , where $p$ is the winding number of $\gamma$ , which we now define. There exists a unique lift $\widehat{\beta}$ of the billiard map $\beta$ to the closure of the universal cover $\mathbb{R}\times[0,\pi]$ which is continuous, $\ell$ periodic in first variable and satisfies $\widehat{\beta}(\widetilde{s},0)=(\widetilde{s},0)$ . Given this normalization, for any point $(\widetilde{s},\widetilde{\sigma})\in\mathbb{R}/\ell\mathbb{Z}\times[0,\pi]$ in the lift of a $q$ periodic orbit of $\beta$ , we see that $\widehat{\beta}^{q}(\widetilde{s},\widetilde{\sigma})=(\widetilde{s}+p\ell,\widetilde{\sigma})$ for some $p\in\mathbb{Z}$ . We define this $p$ to be the winding number of the orbit generated by $\pi(\widetilde{s},\widetilde{\sigma})=(s,\sigma)\in\overline{B^{*}\partial\Omega}$ . Even if a point $(s,\sigma)$ generates an orbit which is not periodic in the full phase space but is such that $\pi_{1}(\beta^{q}(s,\sigma))=s$ for some $q\in\mathbb{Z}$ , we can still define a winding number in this case. Such orbits are called $(p,q)$ geodesic loops, or loops for short. For a given periodic orbit, the winding number is independent of which point in the orbit is chosen, so we sometimes write $\omega(\gamma)=\omega(P^{\prime})$ for any $P^{\prime}\in\{P,\beta(P),\cdots,\beta^{q-1}(P)\}$ . For deeper results and a more thorough introduction to the theory of dynamical billiards, we refer the reader to [Tab95], [Sib04], [Kat05], [Pop94] and [PT11].

3.1. Parametrization and Notation

Following the notation in [MM82], we may translate and rotate $\Omega$ by a rigid motion so that there exists a point $p\in\partial\Omega$ at the point $(0,0)$ with unit tangent vector $(1,0)$ . Let $s$ denote the arclength coordinate along $\partial\Omega$ , measured counterclockwise from $p$ . We then define the arclength parametrization

(4)

\displaystyle x(s)=(x_{1}(s),x_{2}(s))=\int_{0}^{s}(\cos\varphi(t),\sin\varphi(t))dt,

where $\varphi$ is the angle made by the tangent to a base point, $dt$ is arclength on $\partial\Omega$ and $\varphi^{\prime}(t)=\kappa(t)$ is the curvature of $\partial\Omega$ . When dealing with successive points of reflection in the length functional, we use the capital letters $S=(s_{1},\cdots,s_{q})$ to denote the arclength coordinates of the reflection points and $x(S)$ to denote the corresponding boundary points $(x(s_{1}),\cdots,x(s_{q}))$ in $\mathbb{R}^{2}$ , with $x_{i}(S)$ being the $i$ th coordinate. Similarly, when dealing with a periodic orbit $\gamma$ , we denote by $S_{\gamma}$ the arclength coordinates of the reflection points and by $\partial\gamma$ the reflection points in $\mathbb{R}^{2}$ .

3.2. Lengths of Orbits

There are two natural functions for studying the length spectrum. If $\gamma$ is a $q$ -periodic orbit of the billiard map, the length of $\gamma$ is then defined to be

(5)

\displaystyle L_{\gamma}=\sum_{i=1}^{q}|x_{i+1}-x_{i}|,

where $x_{i}$ are the reflection points of $\gamma$ on the boundary, $x_{q+1}=x_{1}$ and $|\cdot|$ is the Euclidean distance function.

Definition 3.1.

For an ordered collection of points $S=(s_{1},\cdots,s_{q})\in[0,\ell)^{q}$ , we define the length functional $\mathcal{L}(S)$ to be $\sum_{i=1}^{q}|x_{i+1}(S)-x_{i}(S)|$ , with $x_{q+1}=x_{1}$ . For any $p,q\in\mathbb{Z}_{>0}$ relatively prime, we define the $(p,q)$ loop function $\ell_{p,q}(s)$ to be the length of the locally unique broken geodesic of rotation number $p/q$ based at $x(s)$ , assuming a well defined branch of it exists (for example, when $p=1$ and $q$ is sufficiently large; see Theorem 3.1 in [Vig23b]).

In the same way that $h(s,s^{\prime})=-|x(s)-x(s^{\prime})|$ is a generating function for $\beta$ , we have that $\mathcal{L}(s_{1},\cdots,s_{q})$ is a generating function for $\beta^{q}$ . The same holds $\ell_{p,q}$ , at least locally in phase space.

Definition 3.2.

The length spectrum of $\Omega$ , denoted by $\text{LSP}(\Omega)$ , is the union of all lengths of periodic orbits.

We will also need the notion of degeneracy for periodic orbits in the Section 6, where we apply the method of stationary phase with the length functional as the phase function of an oscillatory integral.

Definition 3.3.

We say that a periodic orbit $\gamma$ is nondegenerate if $\det(1-P_{\gamma})\neq 0$ , or equivalently if $\partial^{2}\mathcal{L}(S_{\gamma})\neq 0$ (see 6 and 7 below).

3.3. Hessian Invariants

Both generating functions $\mathcal{L}$ and $\ell_{p,q}$ have second order invariants associated to them at a critical point corresponding to a periodic orbit $\gamma$ . The first is the Maslov index, which is a universal affine function of $\text{sgn}\partial^{2}\mathcal{L}\mod 8$ coming from the complex exponential of the signature the Hessian upon performing a stationary phase expansion in 8. From a symplectic geometry perspective, they are intersection numbers of vertical fibers upon iteration by $\beta$ , corresponding to the twist property mentioned above. It is related to the Morse index for the length variational problem with periodic boundary conditions. References can be found in [Bia93] and [WYZ23]. From a microlocal point of view, these indices arise as unitary phases when expanding via stationary phase a Lagrangian distribution associated to the quantization of the billiard map. They come from the integer valued exponents in the unitary transition functions for the Arnold-Keller-Maslov line bundle, sections of which, when tensored with a half density, invariantly define the notion of a principal symbol for Fourier integral operators. Their relavence in the trace formula comes from the fact that the wave propagator is a Fourier integral operator microlocally near nonglancing, transversally reflected rays. The nonpositivity of these complex phases will be important for cancelling contributions to the wave trace of different orbits in our subsequent paper.

The other second order invariant is the determinant of $\partial^{2}\mathcal{L}$ . It is shown in [KT91] (Theorem 3, page 67) that

(6)

\displaystyle\det\partial^{2}\mathcal{L}=(-1)^{q+1}|\det(\text{Id}-P_{\gamma})|\prod_{1}^{q}b_{i},

where $P_{\gamma}$ is the linearized Poincaré map of $\beta$ and

(7)

\displaystyle b_{i}=\frac{\partial^{2}|x_{i}(S)-x_{i+1}(S)|}{\partial s_{i}\partial s_{i+1}}=\frac{\cos\vartheta_{i}\cos\vartheta_{i+1}}{|x_{i}(S)-x_{i+1}(S)|}.

The $b_{i}$ above consist of cosines of angles of reflection for the billiard orbit $\gamma$ at which $\mathcal{L}$ is evaluated and $|x_{i}(S)-x_{i+1}(S)|$ is the length of a link connecting the $i$ th and $i+1$ st reflection points in $\gamma$ . The angles are measured from the inward pointing normal vector at the boundary. The Poincaré map will appear naturally in the symplectic prefactor for our stationary phase expansion of the regularized resolvent trace and the product 6 enters via the multiple reflection expansion with the Hankel functions in Proposition 4.5 below.

Remark 3.4.

The use of $\mathcal{L}:\partial\Omega^{q}\to\mathbb{R}$ as a phase function comes naturally from the layer potential formulas in Section 4. In [MM82], [Vig20] and [HZ22], the $(p,q)$ loop function was used instead. Both are generating functions for the $q$ -th iterate of the billiard map and each presents its own difficulties. One disadvantage of $\mathcal{L}$ is that it has singularities on the diagonal which require regularization. It also has a large number of varaiables to keep track of together with Maslov indices. The loop function is simpler in that it depends only on the base point of an orbit, but is in general multivalued and only local branches exist. Here, we primarily use $\mathcal{L}$ as it is globally well defined and works with potentially smaller $q$ .

4. Balian-Bloch Invariants

Let $\Omega$ be a compact, strictly convex domain having smooth boundary and a nondegenerate periodic orbit $\gamma$ of period $q$ and length $L_{\gamma}$ , which we assume to be isolated in the length spectrum. The connection between the length and Laplace spectrums come from the Poisson relation:

\displaystyle\text{SingSupp}\text{Tr}\,\cos t\sqrt{-\Delta}\subset\{0\}\cup\pm\overline{\text{LSP}(\Omega)}.

The lefthand side is the singular support of the distribution 2, which is a purely spectral quantity, while the righthand side is geometric. In the case of bounded planar domains, this result is due to Anderson and Melrose [AM77]. If $\rho\in\mathcal{S}(\mathbb{R})$ is such that $\hat{\rho}=1$ in a neighborhood of $L_{\gamma}$ and $\text{supp}\hat{\rho}\cap\pm\overline{\text{LSP}(\Omega)}=\{L_{\gamma}\}$ , then the regularized resolvent at frequency $k$ is given by

(8)

\displaystyle\int_{0}^{\infty}e^{ikt}\widehat{\rho}(t)w(t)dt=\frac{1}{i}\text{Tr}\rho\ast kG_{\Omega,D}(k,x,y),

where $G_{\Omega,D}$ is the Green’s kernel, i.e. the Schwartz Kernel of the resolvent $R_{\Omega,D}(k)=(-\Delta_{\Omega,D}-k^{2})^{-1}$ on $\Omega$ with Dirichlet boundary conditions. This identity comes from the formula

\displaystyle R_{\Omega,D}(k)=\frac{i}{k}\int_{0}^{\infty}e^{ikt}\cos t\sqrt{-\Delta_{\Omega,D}}dt.

The Dirichlet resolvent $R_{\Omega,D}$ on $\Omega$ extends meromorphically to $\mathbb{C}$ , with poles at the spectrum of the Laplacian and corresponding residues equal to the orthogonal projectors onto finite dimensional eigenspaces. In view of the Paley-Weiner theorem, the regularized resolvent 8 is in fact an entire function of $k$ . Replacing $k$ by $k+i\tau$ or $k+i\tau\log k$ , the maps $\tau\mapsto e^{-\tau t}\widehat{\rho}(t)w(t)$ and $\tau\mapsto k^{-\tau}\widehat{\rho}(t)w(t)$ give rise to continuous families of compactly supported distributions. Consequently, the regularized resolvent trace has a well defined and smooth limit as $\tau\to 0^{+}$ . When $\tau=0$ , we have the smoothed density of states formula

(9)

\displaystyle\int_{0}^{\infty}e^{ikt}\widehat{\rho}(t)w(t)dt\sim\pi\sum_{\pm}\sum_{n=1}^{\infty}\rho(k\pm\lambda_{n}),

where $\lambda_{n}^{2}$ are the eigenvalues of $-\Delta$ , which provides a connection between the high frequency behavior of $-\Delta$ and the wave trace. From positivity of the spectrum, it is clear that the term in the sum with $\rho(k+\lambda_{n})=O(k^{-\infty})$ is asymptotically negligible. The order of a singularity in the wave trace at a particular length is related to the decay properties of the correspondingly regularized resolvent trace. It follows that the expression 8 has an asymptotic expansion in $k$ of the form

(10)

\displaystyle F_{\gamma}(k)\sum_{j=0}^{\infty}\widetilde{B}_{\gamma,j}k^{-j},

where

(11)

\displaystyle F_{\gamma}(k)=(-1)^{\varepsilon_{B}(\gamma)}\frac{c_{0}L_{\gamma}^{\#}e^{ikL_{\gamma}}e^{i\pi m_{\gamma}/4}}{\sqrt{|\det(I-P_{\gamma})}|}

is called the principal symplectic prefactor. Here, $c_{0}$ is a universal constant, $L_{\gamma}=L$ the length of an isolated nondegenerate periodic orbit $\gamma$ , $L_{\gamma}^{\#}$ the primitive period of $\gamma$ , $m_{\gamma}$ the Maslov index of $\gamma$ , $P_{\gamma}$ the linearized Poincaré map and $\varepsilon_{B}(\gamma)$ depends on boundary conditions (see [Zel09]). The coefficients $\widetilde{B}_{\gamma,j}$ are called the principal Balian-Bloch invariants associated to $\gamma$ , sometimes also called wave trace invariants. The resolvent trace expansion was originally investigated by physicists Balian and Bloch in [BB74], [BB71], [BB72]. The invariants $\widetilde{B}_{j}$ are clearly equivalent to usual wave trace invariants $a_{\gamma,k}$ associated to $w(t)$ near the length spectrum in Theorem 4.2 below. They are independent of the choice of $\rho$ subject to the constraints above, since only the singularity of $w(t)$ at $L_{\gamma}$ affects asymptotics of 10 upon expanding via stationary phase expansion the integral 8.

Our end goal is to show that one can find domains $\Omega$ and lengths of high multiplicity, such that several families of corresponding orbits have Maslov factors of opposite signs. We can then perturb $\partial\Omega$ in such a way that the symplectic prefactors together with $\widetilde{B}_{\gamma,j}$ ’s cancel for arbitrarily large $j$ . If the asymptotic expansion 10 is $O(k^{-\infty})$ , then $w(t)$ is $C^{\infty}$ near $L$ . If the decay is only $O(k^{-m})$ for some $m\geq 2$ , then $w\in C^{m-2}$ locally near $L_{\gamma}$ . One difficulty arises from the nonreality of the coefficients $B_{j}$ , whose complex phases may also depend on the corresponding orbit and could potentially negate whatever complex phase appears from the Maslov index. This is known not to be the case for the leading order invariant $B_{0}$ , but not necessarily for $B_{j}$ with $j\geq 1$ .

We introduce the following notation for simplicity:

Definition 4.1.

Define $\mathcal{D}_{B,\gamma}$ and $B_{\gamma,j}$ by the formula

\displaystyle F_{\gamma}(k)\sum_{j=0}^{\infty}\widetilde{B}_{\gamma,j}k^{-j}=\mathcal{D}_{\gamma}(k)\sum_{j=0}^{\infty}{B}_{\gamma,j}k^{-j},

where

\displaystyle\mathcal{D}_{\gamma}(k)=\frac{c_{0}e^{ikL_{\gamma}}e^{i\pi\text{sgn}\partial^{2}\mathcal{L}/4}}{\sqrt{|\det\partial^{2}\mathcal{L}|}}

and the $B_{\gamma,j}$ are related to the $\widetilde{B}_{\gamma,j}$ as explained above.

The $B_{\gamma,j}$ are called modified Balian Bloch invariants and we will often refer to them simply as the Balian-Bloch invariants when there is no risk of confusion. The reason for introducing $\mathcal{D}_{B,\gamma}$ and $B_{\gamma,j}$ is that they appear more naturally when expanding via stationary phase a microlocal parametrix for the resolvent in Section 6. The corresponding oscillatory integral for the regularized resolvent trace produces terms canonically equivalent to $|\det(\text{Id}-P_{\gamma})|$ and $m_{\gamma}$ in 11, but they appear more naturally in terms of the coordinates $x(S)$ introduced in Section 3.

4.1. Relation to wave invariants

We say that the length $L\in\mathbb{R}$ of a periodic orbit $\gamma$ is simple if up to time reversal ( $t\mapsto-t$ ), $\gamma$ is the unique periodic orbit of length $L$ . Without length spectral simplicity, there is no way to deduce Laplace spectral information from the length spectrum alone. It is shown in [PS92] that generically, smooth convex domains have simple length spectrum and only nondegenerate periodic orbits. In that case, the following theorem holds:

Theorem 4.2 ([GM79], [PS17]).

Assume $\gamma$ is a nondegenerate periodic billiard orbit in a bounded, strictly convex domain with smooth boundary and $\gamma$ has length $L$ which is simple. Then near $L$ , the even wave trace has an asymptotic expansion of the form

\displaystyle\text{Tr}\cos t\sqrt{-\Delta}\sim\operatorname{Re}\left\{a_{\gamma,0}(t-L+i0)^{-1}+\sum_{k=0}^{\infty}a_{\gamma,k}(t-L+i0)^{k}\log(t-L+i0)\right\},

where the coefficients $a_{\gamma k}$ are wave invariants associated to $\gamma$ . The leading order term is given by

\displaystyle a_{\gamma,0}=\frac{e^{i\pi m_{\gamma}/4}L_{\gamma}^{\#}}{|\det(I-P_{\gamma})|^{1/2}}.

Putting the formula from 4.2 into 8, we recover the asymptotics in the Balian-Bloch expansion. An algebraic formula for one set of invariants in terms of the other can be found using basic identities in Fourier analysis.

4.2. Layer Potentials

Here we fix a spectral parameter $\lambda=k+i\tau\in\mathbb{C}_{+}$ with $\tau>0$ and follow [Tay13], [Zel09] in describing potential theory for the resolvent. As before, we let $ds$ denote arclength along the boundary $\partial\Omega$ . From potential theory in the plane (see Proposition 4.7), we have the following formula for the interior Dirichlet resolvent:

(12)

\displaystyle R_{\Omega,D}(\lambda)=\mathds{1}_{\Omega}\left(R_{0}(\lambda)-2\mathcal{D}\ell(\lambda)(\text{Id}+\mathcal{N}(\lambda))^{-1}r_{\partial\Omega}\mathcal{S}\ell^{\dagger}(\lambda)\right)\mathds{1}_{\Omega},

where

	$\displaystyle\mathcal{S}\ell(\lambda)f(x)=\int_{\partial\Omega}G_{0}(\lambda,x,y)f(y)ds(y),\qquad x\in\mathbb{R}^{2}\backslash\partial\Omega$
	$\displaystyle\mathcal{D}\ell(\lambda)f(x)=\int_{\partial\Omega}\partial_{\nu_{s}}G_{0}(\lambda,x,y)f(y)ds(y),\qquad x\in\mathbb{R}^{2}\backslash\partial\Omega$

are the single and double layer potentials at frequency $\lambda$ and the boundary operator $\mathcal{N}(\lambda)$ is given by

(13)

\displaystyle\mathcal{N}(\lambda)f(x)=2\int_{\partial\Omega}\partial_{\nu_{s^{\prime}}}G_{0}(\lambda,x,y)f(y)ds(y).\qquad x\in\partial\Omega

Here, $G_{0}(\lambda,x,y)$ is the free Green’s function, i.e. the Schwartz kernel of the free outgoing resolvent on $\mathbb{R}^{2}$ at frequency $\lambda$ , evaluated at $x\in\mathbb{R}^{2}$ and $y=y(s)\in\partial\Omega$ . The $\pm$ notation indicates limits taken from within $\Omega$ (resp. $\Omega^{c}$ ). The following jump formula holds across $\partial\Omega$ :

\displaystyle\mathcal{D}\ell(\lambda)f_{\pm}=\frac{1}{2}(\pm\text{Id}+\mathcal{N}(\lambda))f.

Remark 4.3.

Note that the signs of the double and boundary layer potentials depend on the choice of unit normal to the boundary. We will use the outward pointing unit normal moving forward, for both the interior and exterior domains.

The free Green’s function on $\mathbb{R}^{2}$ is given by

(14)

\displaystyle G_{0}(\lambda,x,y)=\frac{i}{4}H_{0}^{(1)}(\lambda|x-y|),

where $H_{\nu}^{(1)}$ is a Hankel function of the first kind (of order $\nu$ ) given by $J_{\nu}^{(1)}+iY_{\nu}^{(1)}$ ; $J_{\nu}^{(1)},Y_{\nu}^{(1)}$ are Bessel functions of the first and second kind respectively. The free Greens function solves the problem

\displaystyle(-\Delta_{\mathbb{R}^{2}}-\lambda^{2})G_{0}(\lambda,x,y)=\delta_{0}(x-y).

Remark 4.4.

In even dimensions, the free resolvent $R_{0}(\lambda)$ is only holomorphic on the Riemann surface of the logarithm rather than all of $\mathbb{C}$ , so taking $\lambda=k+i\tau\in\mathbb{C}$ with $\operatorname{Im}\lambda>0$ allows us to consider only the the outgoing resolvent corresponding to $H_{0}(\lambda|x-y|)$ as opposed to $H_{0}(-\lambda|x-y|)$ , since it is bounded on $L^{2}(\mathbb{R}^{2})$ . After using layer potentials to obtain an explicit parametrix for the interior Dirichlet resolvent on $\Omega$ , we can let $\tau\to 0$ using meromorphy of $R_{\Omega,D}$ and holomorphy of the regularized resolvent trace.

Proposition 4.5.

For $\Omega\subset\mathbb{R}^{2}$ with the Euclidean Laplacian, we have

(15)

\displaystyle\begin{split}\mathcal{N}(\lambda,x(s),x(s^{\prime}))\sim\lambda^{1/2}e^{i\lambda|x(s)-x(s^{\prime})|+3\pi i/4}|x(s)-x(s^{\prime})|^{-1/2}\cos\vartheta a_{1}(\lambda|x(s)-x(s^{\prime})|),\end{split}

where $\vartheta$ is the angle made by the link with the interior unit normal at $x(s^{\prime})$ and $a_{1}$ is a semiclassical symbol having the asymptotic expansion

	$\displaystyle a_{1}(z)$	$\displaystyle\sim\sum_{m=0}^{\infty}i^{m}c_{m}z^{-m},\qquad c_{0}=\frac{1}{\sqrt{2\pi}},$
	$\displaystyle c_{m}$	$\displaystyle=\frac{(4-1^{2})(4-3^{2})\cdots(4-(2m-1)^{2})}{m!8^{m}\sqrt{2\pi}}.$

The symbol $a_{1}(\lambda|x(s)-x(s^{\prime})|)$ is holomorphic in $\lambda$ on $\mathbb{C}\backslash(-\infty,0]$ . Tor any fixed $\varrho\in(0,1/2)$ and $\chi\in C_{c}^{\infty}(\mathbb{R})$ supported on $[-1,1]$ , $(1-\chi(k^{\varrho}z))a_{1}((k+i\tau)z)$ belongs to the symbol class $S_{\varrho}^{0}(\mathbb{R})$ :

\displaystyle\left|\partial_{z}^{\alpha}\partial_{k}^{\beta}\left((1-\chi(k^{\varrho}z))a_{0}((k+i\tau)z)\right)\right|\leq C_{\alpha,\beta,\varrho,\tau,K}\langle k\rangle^{\varrho|\alpha|-|\beta|}

for all $K\subset\subset\mathbb{R}$ compact and $\alpha,\beta,\tau$ .

Proof.

Using the well known formula

\displaystyle\frac{d}{dz}H_{\nu}^{(1)}(z)=\frac{\nu H_{\nu}^{(1)}(z)}{z}-H_{\nu+1}^{(1)}(z),

and differentiating $G_{0}$ to obtain the boundary layer operator, we obtain

\displaystyle\mathcal{N}(\lambda,x(s),x(s^{\prime}))=-\frac{i\lambda}{2}\cos\vartheta H_{1}^{(1)}(\lambda|x(s)-x(s^{\prime})|).

The $\cos\vartheta$ term comes from the calculation

\displaystyle\nabla_{y}|x-y|=\frac{y-x}{|x-y|},

which is the unit vector in the direction of the link $\overline{xy}$ . The asymptotics of $H_{1}^{(1)}(\lambda|x(s)-x(s^{\prime})|)$ , including the coefficients $c_{m}$ and the phase $-3\pi/4$ , are given in NIST:

\displaystyle\frac{1}{2}H_{1}^{(1)}(z)\sim\frac{e^{iz-3\pi i/4}}{\sqrt{z}}\sum_{m}i^{m}c_{m}z^{-m}.

Combining with $-i$ gives a total phase of $+3\pi i/4$ . That $a_{1}\in S_{0}^{0}$ follows immediately from the asymptotic expansion in the statement of the proposition. Differentiating term by term gives

\displaystyle\left|\partial_{z}^{\alpha}\partial_{k}^{\beta}a_{0}(kz)\right|\leq\sum_{m=\max{\{\alpha,\beta\}}}^{\infty}c_{m}\frac{(m!)^{2}}{((m-\alpha)!(m-\beta)!)}k^{-m-\beta}z^{-m-\alpha}\lesssim\langle k\rangle^{-\beta}z^{-\alpha},

as $k\to\infty$ , $z$ away from $0$ . Here, $\langle k\rangle=(1+k^{2})^{1/2}$ is the Japanese bracket. Differentiating $(1-\chi(k^{\varrho}z))$ clearly shows that it is in $S_{\varrho}^{0}$ , so $(1-\chi(k^{\varrho}z))a_{1}(kz)\in S_{\varrho}^{0}$ . ∎

Remark 4.6.

The cutoff factor was introduced to localize asymptotics away from zero, where the asymptotic expansion isn’t useful. The introduction of the parameter $\varrho$ is important for the decomposition of $\mathcal{N}$ in Proposition 4.10 below, in which the resolvent is decomposed into a sum of microlocal (homogeneous) pseudodifferential operators and semiclassical Fourier integral operators at scale $\hbar=1/k$ . This decomposition is only valid when working with symbol classes which have $\varrho>0$ .

In calculations below, we will use the function $a_{1}$ in formula 4.5 to avoid extra combinatorial constants while keeping track of complex phases.

4.3. Reduction to the boundary

We now derive the way in which layer potentials translate the interior problem to one defined only on the boundary, following the presentation in [Zel04a]. To do so, we adapt the Poisson integral solution of the boundary value problem considered in [Tay13] to the Dirichlet/Neumann resolvents:

Proposition 4.7.

The Dirichlet resolvent on $\Omega$ and the Neumann resolvent on the exterior domain $\Omega^{c}$ are given by

	$\displaystyle R_{\Omega,D}(\lambda)=$	$\displaystyle\mathds{1}_{\Omega}\left(R_{0}(\lambda)+2\mathcal{D}\ell(\lambda)(\text{Id}+\mathcal{N}_{\Omega}(\lambda))^{-1}r_{\partial\Omega}R_{0}(\lambda)\right)\mathds{1}_{\Omega},$
	$\displaystyle R_{\Omega^{c},N}(\lambda)=$	$\displaystyle\mathds{1}_{\Omega^{c}}\left(R_{0}(\lambda)-2\mathcal{S}\ell(\lambda)(\text{Id}-\mathcal{N}_{\Omega^{c}}^{\#}(\lambda)))^{-1}r_{\partial\Omega^{c}}\partial_{\nu_{\Omega^{c}}}R_{0}(\lambda)\right)\mathds{1}_{\Omega^{c}},$

where

	$\displaystyle\mathcal{N}_{\Omega}^{\#}g(x)=2\int_{\partial\Omega}\partial_{\nu_{x,\Omega}}G_{0}(\lambda,x,y)g(y)ds(y),\qquad x\in\partial\Omega$
	$\displaystyle\mathcal{N}_{\Omega^{c}}^{\#}g(x)=2\int_{\partial\Omega^{c}}\partial_{\nu_{x,\Omega^{c}}}G_{0}(\lambda,x,y)g(y)ds(y).\qquad x\in\partial\Omega^{c}$

Proof.

We begin with solutions for the boundary value problems given in [Tay13], keeping in mind our convention on exterior normals. For brevity, put $P(\lambda)=-\Delta-\lambda^{2}$ :

	$\displaystyle\begin{cases}P(\lambda)u(x)=0,\\ r_{\partial\Omega}u=f,\end{cases}\implies u=-2\mathcal{D}\ell(\lambda)(\text{Id}+\mathcal{N}(\lambda))^{-1}f,\qquad(\text{interior Dirichlet})$
	$\displaystyle\begin{cases}P(\lambda)v(x)=0,\\ r_{\partial\Omega^{c}}\partial_{\nu}v=\varphi,\end{cases}\implies v=2\mathcal{S}\ell(\lambda)(\text{Id}-\mathcal{N}_{\Omega^{c}}^{\#}(\lambda))^{-1}\varphi.\qquad(\text{exterior Neumann})$

Notice that $\mathcal{N}_{\Omega}^{\#}$ is defined with the outward pointing normal for $\Omega^{c}$ , and hence the inward pointing normal for $\Omega$ . To obtain the resolvents, set $\widetilde{G}=R_{\Omega,D}(\lambda)-R_{0}(\lambda)$ and $\widetilde{\widetilde{G}}(\lambda)=R_{\Omega^{c},N}(\lambda)-R_{0}(\lambda)$ . These solve the interior Dirichlet and exterior Neumann problems with $f=-r_{\partial\Omega}R_{0}(\lambda)$ and $\varphi=-r_{\partial\Omega}\partial_{\nu}R_{0}(\lambda)$ . Applying the layer potentials to each and adding back $R_{0}$ completes the proof. ∎

To combine the interior and exterior resolvents into a similar form, note that for the exterior problem, the outer normal points in the opposite direction and $\mathcal{N}_{\Omega}(\lambda)=-\mathcal{N}_{\Omega^{c}}(\lambda)$ . The symmetry $G_{0}(\lambda,x,y)=G_{0}(\lambda,y,x)$ coming from the explicit form of the resolvent in 14 implies that $\mathcal{N}_{\Omega}^{\#}=\mathcal{N}_{\Omega}^{\dagger}$ and $\mathcal{N}_{\Omega^{c}}^{\#}=\mathcal{N}_{\Omega^{c}}^{\dagger}$ . We also have that $(r_{\partial\Omega^{c}}\partial_{\nu}R_{0}(\lambda))^{\dagger}=\mathcal{D}\ell(\lambda)$ and $r_{\partial\Omega}R_{0}(\lambda)=\mathcal{S}\ell^{\dagger}(\lambda)$ . The free resolvent $R_{0}(\lambda)$ is formally self adjoint, so taking the transpose of $R_{\Omega^{c},N}(\lambda)$ and putting in cutoffs to the interior/exterior, we obtain

(16)

\displaystyle\begin{split}R_{\Omega,D}(\lambda)=&\mathds{1}_{\Omega}\left(R_{0}(\lambda)+2\mathcal{D}\ell(\lambda)(\text{Id}+\mathcal{N}(\lambda))^{-1}r_{\partial\Omega}R_{0}(\lambda)\right)\mathds{1}_{\Omega},\\ R_{\Omega^{c},N}^{\dagger}(\lambda)=&\mathds{1}_{\Omega^{c}}\left(R_{0}(\lambda)+2\mathcal{D}\ell(\lambda)(\text{Id}+\mathcal{N}(\lambda))^{-1}r_{\partial\Omega}R_{0}(\lambda)\right)\mathds{1}_{\Omega^{c}},\end{split}

where $r_{\partial\Omega}$ is the operator restricting to the boundary and now all operators are defined using the outward pointing normal from $\Omega$ . Adding the two and subtracting the free resolvent yields

(17)

\displaystyle\text{Tr}\,\left(R_{\Omega,D}(\lambda)+R_{\Omega^{c},N}^{\dagger}(\lambda)-R_{0}(\lambda)\right)=2\text{Tr}\,\left(r_{\partial\Omega}R_{0}(\lambda)(\mathds{1}_{\Omega}^{2}+\mathds{1}_{\Omega^{c}}^{2})\mathcal{D}\ell(\lambda)(\text{Id}+\mathcal{N}(\lambda))^{-1}\right),

where we commuted $r_{\partial\Omega}R_{0}(\lambda)$ through the trace to obtain a kernel on $\partial\Omega\times\partial\Omega$ . Note that the off diagonal terms are trace free. Observe also that the leftmost factor in the righthand side of 17 has the property that

	$\displaystyle 2r_{\partial\Omega}R_{0}(\lambda)\mathcal{D}\ell(\lambda)$	$\displaystyle=2r_{\partial\Omega}R_{0}(\lambda)\circ\partial_{\nu_{2}}R_{0}(\lambda)r_{\partial\Omega}$
		$\displaystyle=2r_{\partial\Omega}\frac{1}{2\lambda}\frac{d}{d\lambda}\partial_{\nu_{2}}R_{0}(\lambda)r_{\partial\Omega}=2\frac{1}{2\lambda}\frac{d}{d\lambda}\left(\frac{1}{2}\mathcal{N}(\lambda)\right),$

where we used the identity

\displaystyle\frac{d}{d\lambda}R_{0}(\lambda)=2\lambda R_{0}^{2}(\lambda),

together with the fact that $\partial_{\nu_{2}}$ is differentiating in only one of the variables. Hence, the righthand side of 17 adds up to

(18)

\displaystyle\frac{1}{2\lambda}\text{Tr}\,\frac{d}{d\lambda}\log(\text{Id}+\mathcal{N}(\lambda))=\frac{1}{2\lambda}\frac{d}{d\lambda}\log\det(\text{Id}+\mathcal{N}(\lambda)).

It is shown in [Zwo98], [Sjö97] and [Chr17] that the regularized exterior resolvent trace $\text{Tr}(R_{\Omega^{c},N}(\lambda)-R_{0}(\lambda))$ admits an asymptotic expansion in negative powers of $\lambda$ for each periodic orbit in $\Omega^{c}$ . By convexity of $\Omega$ , there are only gliding orbits along $\partial\Omega$ and no transversally reflected orbits in the exterior domain. In particular, for $\widehat{\rho}$ supported away from $\mathbb{Z}|\partial\Omega|$ , we have

\displaystyle\int_{0}^{\infty}e^{itk}\widehat{\rho}(t)\cos t\sqrt{-\Delta_{\Omega^{c},N}}=O(k^{-\infty}).

Convolving 18 with the test function $-i\lambda\rho(\lambda)$ gives

(19)

\displaystyle\sum_{\gamma}\mathcal{D}_{\gamma}(k)\sum_{j=0}^{\infty}{B}_{\gamma,j}k^{-j}\sim\frac{1}{2i}\int_{-\infty}^{\infty}\rho(k-\mu)\frac{\partial}{\partial\mu}\log\det(\text{Id}+\mathcal{N}(\mu+i\tau))^{-1}d\mu.

Verification that the relevant functional analytic properties are satisfied to make the above formal computations legitimate (e.g. trace class, Fredholm, etc.) are contained in [Zel04a] and [Zel04c].

4.4. Regularization

As in the introduction, let $\widehat{\rho}\in C_{c}^{\infty}(\mathbb{R})$ be identically equal to $1$ in a neighborhood of an isolated length $L\in\text{LSP}(\Omega)$ , and zero outside of a small interval $(L-\varepsilon,L+\varepsilon)$ such that $\text{supp}\widehat{\rho}\cap\overline{\text{LSP}}=\{L\}$ . Convolving the regularized trace above with $\rho$ and expanding $\mathcal{N}(\lambda)$ in a finite Neumann series yields the following formula connecting the resolvent trace asymptotics to billiard dynamics:

Proposition 4.8 ([Zel09], Proposition 3.6).

Let $\chi_{\partial\gamma}(x,\operatorname{Re}(\lambda)^{-1}D_{x})$ be a semiclassical pseudodifferential operator microlocalizing in phase space near the projection of an orbit $\gamma$ having length $L$ to $B^{*}\partial\Omega$ . Then, for each $J\in\mathbb{N}^{+}$ , there exists $M_{0}(J)$ such that

(20)

\displaystyle\begin{split}\text{Tr}\int_{\mathbb{R}}\rho(k-\zeta)N^{M_{0}+1}(\zeta+i\tau)(\text{Id}+\mathcal{N}(\zeta+i\tau))^{-1}\mathcal{N}^{\prime}(\zeta+i\tau)\chi_{\partial\gamma}(\zeta)d\zeta=O(k^{-(J+1)}).\end{split}

In particular,

(21)

\displaystyle\begin{split}\frac{1}{2i}\sum_{M=0}^{M_{0}}(-1)^{M}\text{Tr}&\int_{\mathbb{R}}\rho(k-\zeta)\mathcal{N}^{M}(\zeta+i\tau)\mathcal{N}^{\prime}(\zeta+i\tau)\chi_{\partial\gamma}(\zeta)d\zeta\\ &=\mathcal{D}_{B,\gamma}(k+i\tau)\sum_{j=0}^{J}B_{\gamma,j}k^{-j}+O(k^{-(J+1)}).\end{split}

This reduces to the computation of $B_{\gamma,j}$ to a calculation of boundary integrals corresponding to powers of $\mathcal{N}$ . The boundary operator $\mathcal{N}$ does not have small operator norm and hence the infinite series isn’t actually convergent. However, the remainder can still be made arbitrarily small for sufficiently large $M$ , which is shown in [Zel04c]. This explains that while not an asymptotic expansion in the usual case, the Neuman series still provides an effective algorithm for computing the coefficients $B_{j}$ . In [BB74], [BB71], [BB72], the physicists Balian and Bloch refer to the sum above as the multiple reflection expansion. The kernel $\mathcal{N}^{M}(\lambda,s,s^{\prime})$ is of the form

	$\displaystyle\mathcal{N}^{M}(\lambda,x(s),x(s^{\prime}))$	$\displaystyle=\int_{\partial\Omega^{M-1}}\mathcal{N}(\lambda,s,s_{1})\mathcal{N}(\lambda,x(s_{1}),x(s_{2}))\cdots\mathcal{N}(\lambda,x(s_{M-1}),x(s^{\prime}))ds_{1}\cdots ds_{M-1}$
		$\displaystyle=\int_{\partial\Omega^{M-1}}e^{ik\mathcal{L}(x(s),x(s_{1}),\cdots,x(s_{M-1}),s^{\prime})}a(\lambda,s,s_{1},\cdots,s_{M-1},s^{\prime})ds_{1}\cdots ds_{M-1},$

for the semiclassical amplitude

(22)

\displaystyle a(\lambda,s,s_{1},\cdots,s_{M-1},s^{\prime})=\lambda^{M/2}e^{3\pi iM/4}\prod_{i=1}^{M}\frac{\cos\vartheta_{i}}{|x(s_{i-1})-x(s_{i})|^{\frac{1}{2}}}a_{1}(\lambda|x(s_{i-1})-x(s_{i})|).

The term multiple reflection refers to the phase function, which is a sum of consecutive generating functions for $\beta$ and hence generates the $M$ -fold iterate of the billiard map. Here, $a_{1}$ is the semiclassical symbol in Proposition 4.5 and we use the convention that $s=s_{0},s^{\prime}=s_{M}$ . Integrating by parts the formula in Proposition 4.8 and reindexing, one obtains

\displaystyle\begin{split}\frac{1}{2i}\sum_{M=1}^{M_{0}+1}\frac{(-1)^{M}}{M}\text{Tr}&\int_{\mathbb{R}}\rho^{\prime}(k-\zeta)\mathcal{N}^{M}(\zeta+i\tau)d\zeta\\ &=\mathcal{D}_{B,\gamma}(k+i\tau)\sum_{j=0}^{J}B_{\gamma,j}k^{-j}+O(k^{-(J+1)}).\end{split}

With the understanding that the following is not a true asymptotic expansion, but rather an algorithm for computing the wave invariants $B_{\gamma,j}$ as above, we write

(23)

\displaystyle\begin{split}\mathcal{D}_{B,\gamma}(k+i\tau)\sum_{j=0}^{\infty}B_{\gamma,j}k^{-j}\sim^{*}\frac{1}{2i}\sum_{M=1}^{\infty}\frac{(-1)^{M}}{M}\text{Tr}\int_{\mathbb{R}}\rho^{\prime}(k-\zeta)\mathcal{N}^{M}(\zeta+i\tau)d\zeta.\end{split}

The notation $\sim^{*}$ is meant to emphasize the fact that the expansion is in the sense of Proposition 4.8. The operator $\mathcal{N}$ is both a microlocal pseudodifferential operater and a semiclassical Fourier integral operator. It has microlocally homogeneous singularities near the diagonal and away from the diagonal, it is a semiclassical Fourier integral operator quantizing the billiard map with semiclassical parameter $\hbar=1/k$ . Hence, some care must be taken to apply the $\Psi$ DO and FIO calculi. In [HZ04], Hassell and Zelditch introduced a microlocal decomposition $\mathcal{N}=\mathcal{N}_{0}+\mathcal{N}_{1}$ ,

	$\displaystyle\mathcal{N}_{0}(x(s),x(s^{\prime}))$	$\displaystyle=\chi(k^{\rho}\|x(s)-x(s^{\prime}))\mathcal{N}(\lambda,x(s),x(s^{\prime})),$
	$\displaystyle\mathcal{N}_{1}(x(s),x(s^{\prime}))$	$\displaystyle=(1-\chi(k^{\rho}\|x(s)-x(s^{\prime})\|))\mathcal{N}(\lambda,x(s),x(s^{\prime})),$

where $\mathcal{N}_{0}\in\Psi^{-1}(\partial\Omega)$ and $\mathcal{N}_{1}\in I_{1/k}^{0}(\partial\Omega\times\partial\Omega;\Lambda_{\beta})$ is a Lagrangian distribution of order zero which is the Schwartz kernel of a semiclassical Fourier integral operator quantizing the billiard map. Here, $\Lambda_{\beta}$ is the twisted graph of the billiard map:

\displaystyle\Lambda_{\beta}=\{(s,s^{\prime},\sigma,-\sigma^{\prime}):(s,\sigma),(s^{\prime},\sigma^{\prime})\in B^{*}(\partial\Omega),\beta(s,\sigma)=(s^{\prime},\sigma^{\prime})\},

thought of as a Lagrangian submanifold of the product cotangent bundle. In this case, one has

\displaystyle\mathcal{N}^{M}=(\mathcal{N}_{0}+\mathcal{N}_{1})^{M}=\sum_{\sigma:\mathbb{Z}_{M}\to\{0,1\}}\mathcal{N}_{\sigma(0)}\mathcal{N}_{\sigma{(1)}}\cdots\mathcal{N}_{\sigma(M)}.

We denote these compositions by $\mathcal{N}_{\sigma}$ and write $|\sigma|=\#\{i:\sigma(i)=0\}$ . The composition rules for Fourier integral operators tell us that the product $\mathcal{N}_{0}\mathcal{N}_{1}$ has canonical relation contained in $\Delta_{\partial\Omega}\circ\Lambda_{\beta}=\Lambda_{\beta}$ , and the fact that $\mathcal{N}_{0}\in\Psi^{-1}$ reduces the order of the symbol by $1$ . Hence, $\mathcal{N}_{\sigma}$ has order $-|\sigma|$ and quantizes the $M-|\sigma|^{th}$ iterate of the billiard map. As we have microlocalized near an orbit $\gamma$ with primitive period $q$ , only the terms with $M-|\sigma|=rq$ for $r\in\mathbb{Z}$ potentially contribute to the leading order asymptotics near the lengths $rL$ of a $p/q$ periodic orbit iterated $r$ times. This is shown rigorously in [Zel09], where it is also demonstrated that for any specified order $R$ and sufficiently large corresponding $M$ , the terms $\mathcal{N}_{0}^{M}$ do not contribute to the regularized resolvent trace asymptotics modulo $O(k^{-R})$ . In addition, it is shown there that if $L=L_{\gamma}$ is simple, then the trace is unchanged modulo $O(k^{-\infty})$ when the semiclassical cutoff in 4.8 is removed. The presence of a microlocal cutoff is what is meant in Theorem 1.1 regarding the contributions to the resolvent trace of a specific orbit $\gamma$ . The reason for a small remainder in Proposition 4.8 is that after microlocalizing near the projection of $\gamma$ to $B^{*}(\partial\Omega)$ and convolving with $\rho$ , $\widehat{\rho}$ being supported near $L_{\gamma}$ , $\mathcal{N}_{\sigma}$ can have at most $q$ factors of $\mathcal{N}_{1}$ without making the phase nonstationary. For large enough $M$ , at least $M-q\gtrsim R$ factors of $\mathcal{N}_{0}$ must enter the composition, each of which reduces the order by $1$ .

The only possible terms contributing to $B_{\gamma,j}$ have $q\leq M\leq q+j$ ; if $M$ were larger, the terms would be at most of order $k^{-j-1}$ and if $M$ were smaller than $q$ , the phase would be nonstationary on the support of $\widehat{\rho}$ , in which case the asymptotics would be $O(k^{-\infty})$ . Moreover, we will see later that only the compositions with $|\sigma|=0$ generate useful terms containing

(1)

highest powers of a yet to be specified deformation parameter $\varepsilon^{-1}$ and
(2)

highest order derivatives of the boundary curvature upon performing a stationary phase expansion.

4.5. Principal symbol calculation

To compute 23, we use the Fourier inversion formula to see that

(24)

\displaystyle\begin{split}&\text{Tr}\,\int\rho^{\prime}(k-\zeta)\mathcal{N}_{\sigma}(\zeta+i\tau)d\zeta=\frac{1}{2\pi}\int_{{\partial\Omega}_{S}^{M}}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}e^{i(k-\zeta)t}it\hat{\rho}(t)\\ &\times\prod_{i=1}^{M}e^{i(\zeta+i\tau)|x(s_{i})-x(s_{i+1})|+3\pi i/4}\chi_{i}(\zeta^{\varrho}|x(s_{i})-x(s_{i+1})|)\\ &\times\sqrt{\frac{(\zeta+i\tau)}{|x(s_{i})-x(s_{i+1})|}}a_{1}((\zeta+i\tau)|x(s_{i})-x(s_{i+1})|)\cos\vartheta_{i}dtd\zeta dS,\end{split}

where $S=(s_{1},\cdots,s_{M})$ , $s_{M+1}=s_{1}$ and

\displaystyle\chi_{i}(z)=\begin{cases}\chi(z)&\sigma(i)=0\\ 1-\chi(z)&\sigma(i)=1,\end{cases}

for a cutoff $\chi$ equal to $1$ on $[-1/2,+1/2]$ and zero outside of $[-1,+1]$ (cf. Proposition 4.5). Changing variables $\zeta\mapsto\zeta/k$ , formula 24 can be written as the oscillatory integral

(25)

\displaystyle\begin{split}k\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{{\partial\Omega}_{S}^{M}}e^{ik\Phi({\zeta}+i{\tau},S,{t})}B_{\sigma}(k{\zeta}+i{\tau},t,S)d{\zeta}d{t}dS,\end{split}

where the phase is given by

(26)

\displaystyle\Phi({\zeta}+i{\tau},S,{t})=(1-\zeta)t+{\zeta}\mathcal{L}(S)

and the amplitude by

(27)

\displaystyle\begin{split}B_{\sigma}(\lambda,S)=&\frac{e^{3M\pi i/4}}{2\pi}it\hat{\rho}({t})\lambda^{\frac{M}{2}}e^{-\tau\mathcal{L}(S)}\prod_{i=1}^{M}\widetilde{\chi}_{i}(\lambda,S)\\ &\times\frac{\cos\vartheta_{i}}{|x(s_{i})-x(s_{i+1})|^{1/2}}a_{1}(\lambda|x(s_{i})-x(s_{i+1})|),\end{split}

where we use the notation

\displaystyle\widetilde{\chi}_{i}(\lambda,S)=\begin{cases}\chi(\operatorname{Re}(\lambda)^{\varrho}|x(s_{i})-x(s_{i+1})|)&\sigma(i)=0\\ 1-\chi(\operatorname{Re}(\lambda)^{\varrho}|x(s_{i})-x(s_{i+1})|)&\sigma(i)=1.\end{cases}

Recall that $\vartheta_{i}$ is the angle between the interior normal to the boundary and the link $(x(s_{i+1})-x(s_{i}))$ . When we later set $|\sigma|=0$ and $M=q$ to obtain leading order asymptotics, all of the $\chi_{i}$ will be one.

Remark 4.9.

Note as above that a priori, $B_{\sigma}\in S_{\varrho}^{\frac{M}{2}}(\partial\Omega^{M})$ . The condition $\varrho\in(0,\frac{1}{2})$ guarantees that one can integrate out the diagonal singularities corresponding to terms with $\sigma(i)=0$ to obtain a Fourier integral kernel of lower order quantizing $\beta^{M-|\sigma|}$ , as was done in [Zel09]. We will make use of this in Section 6.

An application of stationary phase in the variables ${\zeta},{t}$ of 25 gives the following:

Proposition 4.10 ([Zel09]).

For each $\sigma:\mathbb{Z}_{M}\to\{0,1\}$ , $\text{Tr}\rho^{\prime}\ast\mathcal{N}_{\sigma}(k+i\tau)$ is given modulo $k^{-\infty}$ by

(28)

\displaystyle e^{3M\pi i/4}i\int_{{\partial\Omega}_{S}^{M}}e^{i(k+i\tau)\mathcal{L}(S)}a_{0}^{\sigma}(k+i\tau,S)dS

where

\displaystyle a_{0}^{\sigma}(k+i\tau,S)=\hat{\rho}(\mathcal{L}(S))\left(\mathcal{L}(S)A_{\sigma}(k+i\tau,S)-\frac{1}{i}\frac{\partial}{\partial k}A_{\sigma}(k+i\tau,S)\right)\in S_{\varrho}^{\frac{M}{2}}(\partial\Omega^{M})

and $A_{\sigma}(k+i\tau,S)$ is

\displaystyle(k+i\tau)^{M/2}\prod_{i=1}^{M}\chi_{i}(k+i\tau,S)\frac{\cos\vartheta_{i}}{|x(s_{i})-x(s_{i+1})|^{1/2}}a_{1}((k+i\tau)|x(s_{i})-x(s_{i+1})|).

The complex phase is given by $(3M+2)\pi/4$ and $a_{1}$ is the Hankel amplitude from Proposition 4.5.

Proof.

See Equation 32 in Section 6 below for the stationary phase formula. In our case, the stationary points occur at $t=\mathcal{L},{\zeta}=1$ . The Hessian is given by

\displaystyle\partial_{{\zeta},{t}}^{2}\Phi=\begin{pmatrix}0&-1\\ -1&0\end{pmatrix},

from which we obtain $\langle\partial_{{\zeta},{t}}^{2}\Phi^{-1}\partial_{{\zeta},{t}},\partial_{{\zeta},{t}}\rangle=-2\partial_{{\zeta}}\partial_{{t}}$ . The signature of the Hessian is $0$ and the modulus of its determinant is one. The total prefactor is then $\frac{2\pi}{k}e^{ik\mathcal{L}}$ , the coefficient of which cancels with that from Fourier inversion and the change of variables above. Since the phase is quadratic, only the amplitude is differentiated ( $g=0$ , $\mu=0$ in the notation of formula 32 below). Since the amplitude factors into a locally linear function of $t$ and one of $\zeta$ , there are only two nonzero terms. Near the critical point corresponding to the length $L_{\gamma}$ of a periodic orbit, $\hat{\rho}$ is constant and hence $\partial_{t}$ only lands on $t$ . The $k$ coming from $\partial_{{\zeta}}$ is canceled by the $k^{-1}$ in stationary phase as is the factor of $2\pi$ . The factor $e^{-\tau\mathcal{L}(S)}$ is unchanged by the Hessian operator and can be reinserted into the phase function. ∎

In the next section, we will prove that the $|\sigma|\geq 1$ terms do not contribute to the leading order deformation asymptotics or highest order derivatives in the curvature jet of the boundary. Hence, we can set $\chi_{i}=1$ in the expression for $A_{\sigma}$ in 4.10. As $\mathcal{N}_{0}\in\Psi^{-1}(\partial\Omega)$ , the symbol $a_{0}^{\sigma}$ in Proposition 4.10 can be integrated to yield a symbol in $S_{\varrho}^{\frac{M-3|\sigma|}{2}}(\partial\Omega^{M-|\sigma|})$ ; see Theorem 3.8 in [Zel09]. Moreover, since $a_{1}\in S_{\varrho}^{\frac{M}{2}}(\partial\Omega^{M})$ , $A_{\sigma}^{\prime}\in S_{\varrho}^{\frac{M}{2}-1}$ . Hence, for $|\sigma|=0$ , $M=q$ and $S$ near $S_{\gamma}$ , we have

(29)

\displaystyle a_{0}^{0}(\lambda,S)=\left(\frac{\lambda}{2\pi}\right)^{q/2}\mathcal{L}(S)\prod_{i=1}^{q}\frac{\cos\vartheta_{i}}{|x(s_{i})-x(s_{i+1})|^{1/2}}+O_{S_{\varrho}^{0}}\left(\operatorname{Re}(\lambda)^{\frac{q}{2}-1}\right),

where we used the leading order asymptotics for $a_{1}$ in formula 4.5.

Remark 4.11.

Recall that by formula 6, the product in the leading order expansion of $a_{0}$ is precisely

\displaystyle\prod_{i=1}^{q}\frac{\cos\vartheta_{i}}{|x(s_{i})-x(s_{i+1})|^{1/2}}=\sqrt{\frac{\det\partial^{2}\mathcal{L}}{|\det(1-P_{\gamma})|}},

with $P_{\gamma}$ being the Poincaré map associated to $\gamma$ .

Combining Proposition 4.10 and 23, we obtain the following:

Corollary 4.12.

The coefficients in the regularized resolvent trace expansion can be determined from the oscillatory integrals in Proposition 4.10:

\displaystyle\mathcal{D}_{\gamma}(\lambda)\sum_{j=0}^{\infty}B_{\gamma,j}\operatorname{Re}(\lambda)^{-j}\sim^{*}\frac{1}{2}\sum_{M=1}^{\infty}\sum_{\sigma:\mathbb{Z}_{M}\to\{0,1\}}\frac{e^{-M\pi i/4}}{M}\int_{{\partial\Omega}_{S}^{M}}e^{i\lambda\mathcal{L}(S)}a_{0}^{\sigma}(\lambda,S)dS.

5. Perturbation Thoery

Starting with a domain $\Omega_{0}$ and a degenerate periodic orbit $\gamma$ , we introduce a one parameter family of domains $\Omega_{\varepsilon}$ specified by a smooth function $\mu$ :

\displaystyle\Omega_{\varepsilon}=\{x+\mu(x,\varepsilon)\nu_{x}:x\in\partial\Omega_{0}\},

where $\nu_{x}$ is the outward pointing normal to $\Omega_{0}$ at $x$ . We impose the following constraints on $\Omega_{\varepsilon}$ :

(1)

$\Omega_{0}=\Omega$ .
(2)

$\Omega_{\varepsilon}$ and $\Omega_{0}$ make first order contact at the reflection points $x_{i}\in\partial\gamma$ .
(3)

Near each $x_{i}$ , $\Omega_{\varepsilon}$ is given by $x+\mu(x,\varepsilon)\nu(x)$ , where $x\in\partial\Omega$ and $\nu(x)$ is the unit normal at $x$ .
(4)

$\mu$ is $C^{\infty}$ in all parameters.
(5)

$\det\partial^{2}\mathcal{L}\sim c_{\gamma}\varepsilon$ for $c_{\gamma}\neq 0$ .
(6)

$\partial^{2}\mathcal{L}$ has rank $q-1$ on $\Omega_{0}$ .
(7)

The third and higher order derivatives of $\mathcal{L}$ are $O(1)$ near $\partial\gamma$ .

In particular, the perturbation preserves periodicity of $\gamma$ and removes its degeneracy. It will be important for us that the Hessian of $\mathcal{L}$ on $\partial\Omega_{0}$ be minimally degenerate in the sense that it has rank $q-1$ . Theorem 1.1 applies to any family of domains which satisfy the conditions above. To demonstrate the robustness of our main theorem, we now give a large class of examples to which it applies.

Definition 5.1.

Given a domain $\Omega$ and a $(p,q)$ -periodic orbit $\gamma$ , we say that $\gamma$ satisfies the injectivity condition if its relfection points are all distinct.

Proposition 5.2.

Let $\mathcal{C}_{1}$ be the class of smooth, bounded strictly convex planar domains $\Omega$ which have a degenerate periodic orbit $\gamma$ satisfying the injectivity condition. Within $\mathcal{C}_{1}$ , there is an open and dense set in the $C^{\infty}$ topology on which the Hessian of $\mathcal{L}$ has rank $q-1$ .

Proposition 5.3.

Given a domain $\Omega_{0}$ and a degenerate orbit $\gamma$ having $\text{rank}(\partial^{2}\mathcal{L}(S_{\gamma}))=q-1$ , consider the class $\mathcal{C}_{2}$ of deformations $\mu:\partial\Omega\times[0,\varepsilon_{0})\to\mathbb{R}$ which preserve $\gamma$ and the angles of reflection. Within $\mathcal{C}_{2}$ , those which make $\partial^{2}\mathcal{L}$ nondegenerate at $S_{\gamma}$ in arclength coordinates form an open set. If in addition, the orbit satisfies the injectivity condition, then they are also dense in $\mathcal{C}_{2}$ .

The preceding two propositions follow from the multijet transversality theorem, which is a generalization due to Mather of Thom’s transversality theorem. It is easy to check that for ellipses, the Hessian has rank $q-1$ . However, degenerate homoclinic orbits such as those constructed in [Cal22] could have lower rank Hessians. In [KZ18], it is shown that domains with a rational caustic of rotation number $1/q$ are quantitatively dense in the space of all convex domains. For analytic domains, the density is exponential in $q$ , whereas for domains of finite smoothness, the density is polynomial. For such domains, it is easy to see that the length functional evaluated at reflection points of an orbit tangent to a corresponding rational caustic has Hessian of rank $q-1$ . Together with the cancellations in our forthcoming paper, we see that the any finite order singular support of the wave trace is generically distinct from the length spectrum in the sense above.

5.1. Nonsingular perturbations

We first show that the jet of the length functional is algebraically equivalent to that of the boundary curvature at reflection points and in particular, demonstrate that arbitrarily small perturbations of curvature can result in nondegeneracy. In what follows, write $\det\mathcal{L}_{\varepsilon}^{-1}=M_{\varepsilon}$ .

Proposition 5.4.

If the rank of $\partial^{2}\mathcal{L}$ is $q-1$ , then the condition $\det\partial^{2}\mathcal{L}_{\varepsilon}\neq 0$ imposes a codimension $1$ constraint on the $q$ -fold $0$ -jet of the boundary curvature at reflection points. For each point $x_{i}\in\partial\gamma$ and each $1\leq i\leq q$ , the data $\vartheta_{i}$ , $|x_{i}-x_{i+1}|$ and $\partial_{i_{1},\cdots,i_{m}}^{m}\mathcal{L}$ , with $|i_{j}-i_{k}|\leq 1$ , uniquely determine the $q$ -fold $m-2$ -jet of the curvature at $x_{i}\in\partial\gamma$ .

Proof.

It is shown in [KT91] that

(30)

\displaystyle\begin{split}\delta_{i}&:=\frac{\partial^{2}\mathcal{L}}{\partial^{2}s_{i}}=2\cos\vartheta_{i}\left(\frac{\cos\vartheta_{i}}{\ell_{i}}-\kappa_{i}\right),\\ \alpha_{i}&:=\frac{\partial^{2}\mathcal{L}}{\partial s_{i}\partial s_{i+1}}=\frac{\cos\vartheta_{i}\cos\vartheta_{i+1}}{\ell_{i}},\end{split}

where $\ell_{i}=|x_{i}(S)-x_{i+1}(S)|$ , $\vartheta_{i}$ is the angle of incidence with the inward pointing normal at $x_{i}(S)$ , and $\kappa_{i}$ is the curvature of $\partial\Omega$ at $x_{i}(S)$ . It follows that $\partial^{2}\mathcal{L}$ is tridiagonal and cyclic with $\delta_{i}$ s on the diagonal and $\alpha_{i}$ s on the off diagonals. The Hessian matrix is of the form

(31)

\displaystyle\partial^{2}\mathcal{L}=\begin{pmatrix}\delta_{1}&\alpha_{1}&0&0&\cdots&\alpha_{q}\\ \alpha_{1}&\delta_{2}&\alpha_{2}&0&\cdots&0\\ 0&\alpha_{2}&\delta_{3}&\alpha_{3}&\cdots&0\\ 0&0&\alpha_{3}&\delta_{4}&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ \alpha_{q}&0&0&0&\alpha_{q-1}&\delta_{q}\end{pmatrix}

and the equivalence of curvature and length functional jets follows immediately. Note that the angles of incidence and lengths of links are preserved under deformations which make first order contact with $\partial\Omega_{0}$ at $\partial\gamma$ . However, the curvature $\kappa_{i}$ in the formula for $\delta_{i}$ can change. If the perturbed curvature has an asymptotic expansion of the form

\displaystyle\kappa_{i}(\varepsilon)\sim\sum_{j=0}^{\infty}\mu_{j}^{i}\varepsilon^{j}

for some smooth functions $\mu_{j}^{i}(x)$ , then we have

	$\displaystyle\det\partial^{2}\mathcal{L}_{\varepsilon}$	$\displaystyle=0+\varepsilon\text{Tr}\left(\text{adj}\left(\partial^{2}\mathcal{L}_{0}\right)\partial_{\varepsilon}\partial^{2}\mathcal{L}_{\varepsilon}\big{\|}_{\varepsilon=0}\right)+O(\varepsilon^{2})$
		$\displaystyle=-\text{Tr}(M\mu)\varepsilon+O(\varepsilon^{2}),$

where $\mu$ is the diagonal matrix given by $\mu_{1}^{i}\delta_{ij}$ and $M$ is the adjugate of $\partial^{2}\mathcal{L}_{0}$ consisting of principal minors. For any matrix, its rank is given by the maximal size for which there is a nonzero minor. If $\partial^{2}\mathcal{L}$ has rank $q-1$ , then the rank of $M$ is $1$ . Hence, the condition $\text{Tr}M\mu\neq 0$ imposes a codimension one, linear constraint on the perturbation matrix $\mu$ or equivalently, on the curvature at reflection points. ∎

From the formulas 30 and 31 above, one sees that the rank of $\partial^{2}\mathcal{L}$ is always at least $q-2$ due to the off diagonal terms being nonzero for nonglancing orbits. If it is in fact equal to $q-2$ , then the linearized Poincaré map is equal to the identity. If the rank of $\partial^{2}\mathcal{L}$ is $q-1$ , the linearized Poincaré map has $1$ as an eigenvalue but differs from the identity and if the rank is $q$ , then the orbit is nondegenerate.

5.2. Analysis of the Hessian

Proposition 5.5.

If $\partial^{2}\mathcal{L}$ has rank $q-1$ in arclength coordinates at a degenerate periodic trajectory $\gamma$ in $\Omega_{0}$ , then its adjugate matrix has rank $1$ and its elements are of the form $h_{i_{1},i_{2}}(S,0)=\pm h_{i_{1}}(S,0)h_{i_{2}}(S,0)$ for some smooth functions $h_{i}(S,0)$ . If a deformation $\Omega_{\varepsilon}$ which preserves the points and angles of reflection makes $|\det\partial^{2}\mathcal{L}|\sim c_{\gamma}\varepsilon$ , $c_{\gamma}>0$ , then the $(i_{1},i_{2})$ entry of $\partial^{2}\mathcal{L}^{-1}$ is $\widetilde{h}_{i_{1}i_{2}}=\sigma\cdot\widetilde{h}_{i_{1}}\widetilde{h}_{i_{2}}+O(1)$ with each $\widetilde{h}_{i}$ having the form $(c_{\gamma}\varepsilon)^{-1/2}h_{i}$ for some smooth functions $h_{i}\in C^{\infty}(\mathbb{R}/\ell\mathbb{Z}\times[0,\varepsilon_{0}))$ which are uniformly bounded in $\varepsilon$ . Here, $\sigma$ is $+1$ if closest eigenvalue to zero of $\partial^{2}\mathcal{L}$ is positive and $-1$ otherwise. If $\Omega_{0}$ is a circle, then $h_{i_{1}}=h_{i_{2}}$ for all $1\leq i_{1},i_{2}\leq q$ .

Proof.

From Proposition 5.4 above, we can choose a perturbation $\mu$ such that $\det\partial^{2}\mathcal{L}\sim c_{\gamma}\varepsilon$ , $c_{\gamma}\neq 0$ . Denote by $H=H_{\varepsilon}=\partial^{2}\mathcal{L}$ the Hessian of $\mathcal{L}$ on $\partial\Omega_{\varepsilon}$ . The adjugate matrix $M_{\varepsilon}$ of $\partial^{2}\mathcal{L}$ in $\Omega_{\varepsilon}$ consists of complimentary minors depending on $\|\kappa\|_{C^{0}(\partial\Omega_{\varepsilon})}$ . If the perturbation $\mu(\varepsilon,\cdot)$ makes $\det\mathcal{L}_{\varepsilon}\sim c_{\gamma}\varepsilon\neq 0$ , then the inverse Hessians satisfy

\displaystyle\partial^{2}\mathcal{L}^{-1}\sim\frac{1}{c_{\gamma}\varepsilon}M_{\varepsilon},

where the cofactor matrices are $C^{\infty}(\mathbb{R}^{q^{2}}\times[0,\varepsilon_{0}])$ and uniformly bounded in $\varepsilon$ all the way down to $\varepsilon=0$ . In particular, $\partial^{2}\mathcal{L}^{-1}\sim(c_{\gamma}\varepsilon)^{-1}M_{0}+O(1)$ . Under the rank $q-1$ assumption, $M_{0}$ is rank $1$ , which together with symmetry about the diagonal implies that $\partial^{2}\mathcal{L}^{-1}\sim\sigma\cdot(c_{\gamma}\varepsilon)^{-1}{h}_{i_{1}}{h}_{i_{2}}$ as in the statement of the theorem.

If $\Omega_{0}$ is a circle, then we have action angle variables in which the angles are just the usual arclength coordinates on $\partial\Omega_{0}$ and the actions are the radii of concentric circles to which corresponding orbits are tangent. It is easy to see in this case that $\partial^{2}\mathcal{L}_{0}$ is negative semi-definite with rank $q-1$ and has kernel $\mathbb{R}v$ , for $v=(1,\cdots,1)$ corresponding to rotation along a concentric circle. We claim that each row of $M_{0}$ is a multiple of $v$ . Let $M^{i}$ denote the $i$ th row of $M_{0}$ and choose any $w$ orthogonal to $v$ . Then, then using the definition of $M$ in terms of complimentary minors, we have

\displaystyle\langle M^{1},w\rangle=\det\begin{pmatrix}w_{1}&w_{2}&\cdots&w_{q}\\ H_{0}^{21}&H_{0}^{22}&\cdots&H_{0}^{2q}\\ \vdots&\ldots&\ddots&\vdots\\ H_{0}^{q1}&H_{0}^{q2}&\ldots&H_{0}^{qq}\end{pmatrix},

where $H_{0}^{ij}$ are the entries of $\partial^{2}\mathcal{L}_{0}$ . Since $w\perp\ker H_{0}=\ker H_{0}^{*}$ , $w\in\text{Span}\{H_{0}^{1},H_{0}^{2},\cdots,H_{0}^{q}\}=\text{Span}\{H_{0}^{2},\cdots,H_{0}^{q}\}$ , given that the rows sum to zero and are hence linearly dependent. Therefore, $\langle M^{1},w\rangle=0$ . We have shown that $\text{Span}(v)^{\perp}\subset\text{Span}(M^{1})^{\perp}$ which implies $\text{Span}(M^{1})=\text{Span}(v)$ since they have equal dimensions; the same holds if one replaces $1$ by any $1\leq i\leq q$ . Each row of $M$ is a multiple of $v$ and by symmetry about the diagonal, the multiples must be the same. In particular,

\displaystyle\partial^{2}\mathcal{L}^{-1}\sim\frac{\sigma}{c_{\gamma}\varepsilon}h(1,\cdots,1)^{T}(1,\cdots,1)+O(1),

for a single smooth function $h$ on $\Omega_{0}$ . Since $\partial^{2}\mathcal{L}$ has $q-1$ negative eigenvalues on the disk, $\sigma$ only depends on ellipticity or hyperbolicity of the corresponding perturbed orbit. ∎

6. Stationary Phase and Feynman Diagrams

We now apply the method of stationary phase to the integral appearing in Theorem 4.10, but on the perturbed domain $\Omega_{\varepsilon}$ from Section 5, so as to make the phase nondegenerate. In general, if $\Phi$ has a unique stationary point at $x_{0}$ with nondegenerate Hessian and $u$ is a compactly supported smooth function, the stationary phase formula in [Hör03] reads:

(32)

\displaystyle\int_{\mathbb{R}^{n}}u(x)e^{ik\Phi(x)}dx\sim(2\pi/k)^{n/2}\frac{e^{ik\Phi(x_{0})}e^{i\pi\text{sgn}\partial^{2}\Phi(x_{0})/4}}{|\det\partial^{2}\Phi(x_{0})|^{1/2}}\sum_{j=0}^{\infty}k^{-j}L_{j}u(x_{0}),

where the $L_{j}$ are differential operators of order $2j$ having the form

(33)

\displaystyle L_{j}u(x_{0})=\sum_{\nu-\mu=j}\sum_{2\nu\geq 3\mu}i^{-j}2^{-\nu}\langle\partial^{2}\Phi(x_{0})^{-1}\partial,\partial\rangle^{\nu}(g^{\mu}u(x_{0}))/\mu!\nu!.

Here, $g=\Phi(x)-\Phi(x_{0})-\Phi^{\prime}(x_{0})(x-x_{0})-\Phi^{\prime\prime}(x_{0})(x-x_{0})^{2}$ is the higher order part of the Taylor expansion of $\Phi$ , having polynomial order $\geq 3$ . We in fact need a version of stationary phase which allows for dependence on auxilary parameters. If the phase $\Phi$ and amplitude $u$ depend smoothly on a parameter $\tau\in T$ with $u$ being supported in a compact set which is independent of $\tau$ , then expansion 32 is uniform in $\tau$ (see [Dui96]). There are several asymptotics to keep track of simultaneously: powers of $k$ and $\varepsilon$ together with the number of boundary curvature derivatives.

We use Feynman diagrams to organize the terms hierarchically, inspired by but using a different structure than those in [Zel09].

6.1. Feynman Calculus

We review the procedure first used by Feynman to evaluate asymptotic integrals. We will follow the conventions on notation as in [Axe97].

Definition 6.1.

A graph $\mathcal{G}$ with $V$ vertices and $I$ edges is called a Feynman diagram of order $j$ if

(1)

$I-V=j$ ,
(2)

there are $V$ closed (indistinguishable) vertices, each of which has valency at least $3$ ,
(3)

there is one open, marked vertex with any possible valency, including $0$ (corresponding to an isolated vertex), and
(4)

each edge is labeled by a function $\ell$ , which assigns two indices $(i_{1},i_{2})$ corresponding two its two endpoints.

Remark 6.2.

The quantity $I-V=j$ is closely related to the Euler characteristic $\chi(\mathcal{G})$ , but does not include faces; not all diagrams are assumed to be planar graphs. We will sometimes denote $\mathcal{G}$ by $\mathcal{G}_{\ell}$ to emphasize $\ell$ , the edge labeling.

Definition 6.3.

Let $\mathcal{G}$ be a Feynman diagram of order $j$ with labeling $\ell$ .

•

For the open vertex, we associate the quantity

$\frac{\partial^{v}u}{\partial x_{i_{1}}\cdots\partial x_{i_{v}}}(x_{0}),$

where $i_{1},\cdots,i_{v}$ are the indices labeling the incident edges ( $v$ is the valency of the vertex).
•

To each closed vertex corresponds a factor of

$ik\frac{\partial^{v}\Phi}{\partial x_{i_{1}}\cdots\partial x_{i_{v}}}(x_{0}),$

where again $i_{1},\cdots,i_{v}$ are the indices of the incident edges.
•

For each edge with indices $i_{1},i_{2}$ , we assign the inverse Hessian element

$\frac{i}{k}(\partial^{2}\Phi)_{i_{1}i_{2}}^{-1}:=\frac{i}{k}\widetilde{h}_{i_{1}i_{2}}.$

The Feynman amplitude $F_{\mathcal{G},\ell}$ corresponding to $\mathcal{G}$ and label $\ell$ is given by the product of each factor above.

Proposition 6.4 ([Axe97]).

Let $Z_{k}$ denote the integral 32. As $k\to\infty$ ,

(34)

\displaystyle Z_{k}\sim\left(\frac{2\pi}{k}\right)^{n/2}\frac{e^{ik\Phi(x_{0})}e^{i\pi\text{sgn}\partial^{2}\Phi}}{|\det\partial^{2}\Phi|^{1/2}}\sum_{I=0}^{\infty}\sum_{V=0}^{\infty}\sum_{\mathcal{G},\ell}\frac{F_{\mathcal{G}}}{w(\mathcal{G})},

where the inner sum is over all labeled labeled Feynman diagrams with $V$ internal vertices, $1$ external vertex, and $I$ edges. $\omega(\mathcal{G})$ is the order of the automorphism group of $\mathcal{G}$ , $F_{\mathcal{G},\ell}$ is the Feynman amplitude and $\ell$ assigns to each endpoint of each edge an index in $\{1,\cdots,n\}$ .

We refer to [Axe97] for an elementary proof of the equivalence of 32 and Proposition 6.4. It remains evaluate each of these Feynman amplitudes in coordinates to determine their contributions to $B_{\gamma,j}$ .

6.2. Contributing Graphs

Here, we examine the Feynman diagrams which contribute maximal inverse powers of $\varepsilon$ to the trace expansion in Proposition 4.10. It is clear that the amplitudes $a_{0}^{\sigma}(k+i\tau,S;\varepsilon)$ are uniformly bounded in $\varepsilon$ . The main contributions come from the closed vertices (derivatives of the phase) and their incident edges (inverse Hessians). To keep track of leading order terms in $\varepsilon$ , we must isolate those Feynman diagrams which produce the most inverse Hessians, each of which produces a power of $\varepsilon^{-1}$ .

Lemma 6.5.

For each $j$ , the Feynman amplitudes having maximal products of inverse Hessian elements correspond to the class of $3$ -regular graphs on $2j$ closed vertices with one additional open vertex having valency zero.

Proof.

In the usual non-diagrammatic stationary phase notation, maximal inverse Hessians arise from the terms where $\nu-\mu=j$ and $2\nu=3\mu$ , which yields $\nu=2j$ and $\mu=3j$ . In Proposition 6.4, $V-I=j$ and $\mathcal{G}$ has $I=3j$ edges. For a general graph,

I=\frac{u_{0}}{2}+\frac{1}{2}\sum_{i=1}^{V}v_{i},

where $u_{0}$ is the valency of the open vertex and $v_{i}$ are the valencies of unmarked vertices. As each $v_{i}\geq 3$ by assumption, we have

	$\displaystyle I-\frac{3}{2}V\geq$	$\displaystyle 0,$
	$\displaystyle-\frac{3}{2}I+\frac{3}{2}V=$	$\displaystyle-\frac{3}{2}j,$

from which it follows that $I\leq 3j$ , which is the maximal number of inverse Hessians. This is acheived when $V=2j$ and $I=3j$ , corresponding to each $v_{i}=3$ . If the open vertex were not isolated, consider the graph $\mathring{\mathcal{G}}=\mathcal{G}\backslash\{\text{open vertex and all its incident edges}\}$ . The formula above would then force there to be vertices with valency strictly less than $3$ . ∎

Figure 1. Some admissible Feynman diagrams given by

3

regular graphs on

2j

vertices.

Corollary 6.6.

Modulo an error bounded uniformly in $\varepsilon$ , the Feynman amplitudes for contributing graphs with maximal Hessians consist of $a_{0}(S_{\gamma})(\pm i)^{j}k^{-j}$ and products of terms of the form

\displaystyle(c_{\gamma}\varepsilon)^{-\frac{3}{2}}h_{i_{1}}h_{i_{2}}h_{i_{3}}\sum_{i_{1},i_{2},i_{3}}\partial_{i_{1}i_{2}i_{3}}^{3}\mathcal{L},

for some indices $1\leq i_{1},i_{2},i_{3}\leq q$ and $a_{0}$ the amplitude defined in Proposition 4.10.

6.3. Critical points on the diagonal

The integral in Theorem 4.10 has two types of critical points. There are those on the large diagonals $x_{i}=x_{i+1}$ and regular critical points which are off diagonal.

Lemma 6.7.

[[Zel09]] The amplitude $A_{\sigma}$ of $\mathcal{N}_{\sigma}$ is a semiclassical amplitude of order $-|\sigma|$ .

Hence, the critical points on the large diagonal have lower order and we can focus on the portion of $N$ which quantizes the $M$ -fold billiard map on $S_{+}^{*}\partial\Omega$ .

Corollary 6.8.

The coefficient of $k^{j}$ in the asymptotic expansion of the integral coming from $\mathcal{N}_{\sigma}$ in Theorem 4.10 is $O(\varepsilon^{-3j+|\sigma|})$ .

Proof.

By Lemma 6.7, $\text{Tr}\rho^{\prime}\ast\mathcal{N}_{\sigma}(k)$ has an asymptotic expansion in negative powers of $k$ begining with $k^{-|\sigma|}$ . Hence, the coefficient of $k^{-j}$ corresponds to the $j-|\sigma|$ term in the usual stationary phase expansion 32 and the maximal contribution of inverse Hessians to such a term is $3(j-|\sigma|)\leq 3j$ . ∎

When $|\sigma|>0$ , these terms are lower order in $1/\varepsilon$ and can be included in the remainder so it suffices to only consider the contributions of regular critical points in the limit $\varepsilon\to 0$ . Notice also that the $A_{\sigma}^{\prime}$ term in the amplitude $a_{0}$ appearing in Theorem 4.10 is a semiclassical symbol of order $-1-|\sigma|$ and can be discarded from the leading order asymptotics in $k$ by the same reasoning as in Corollary 6.8.

6.4. Proof of Theorem 1.1

We can now complete the derivation of the regularized resolvent trace formula appearing in Theorem 1.1.

Proof.

We will analyze the expression

\displaystyle\mathcal{D}_{\gamma}(\lambda)\sum_{j=0}^{\infty}B_{\gamma,j}\operatorname{Re}(\lambda)^{-j}\sim^{*}\frac{1}{2}\sum_{M=1}^{\infty}\sum_{\sigma:\mathbb{Z}_{M}\to\{0,1\}}\frac{e^{-M\pi i/4}}{M}\int_{{\partial\Omega}_{S}^{M}}e^{i\lambda\mathcal{L}(S)}a_{0}^{\sigma}(\lambda,S)dS

in Proposition 4.10 by first expanding the amplitude $a_{0}^{\sigma}$ in inverse powers of $k$ from the representation 4.5 and then replacing each term by its stationary phase expansion. Recall that

\displaystyle a_{0}^{\sigma}(k+i\tau,S)=\hat{\rho}(\mathcal{L}(S))\left(\mathcal{L}(S)A_{\sigma}(k+i\tau,S)-\frac{1}{i}\frac{\partial}{\partial k}A_{\sigma}(k+i\tau,S)\right),

with $A_{\sigma}(k+i\tau,S)$ being

\displaystyle(k+i\tau)^{M/2}\prod_{i=1}^{M}\frac{\chi_{i}^{\sigma}(k+i\tau,S)}{|x(s_{i})-x(s_{i+1})|^{1/2}}a_{1}((k+i\tau)|x(s_{i})-x(s_{i+1})|)\cos\vartheta_{i},

and $a_{1}$ the Hankel amplitude which has the asymptotic expansion 4.5. It is shown in [Zel09] (Corollary 3.9 and the proof of Proposition 3.10, together with Lemma 5.8) that

\displaystyle\int_{\partial\Omega^{M}}e^{i(k+i\tau)\mathcal{L}(S)}a_{0}^{\sigma}((k+i\tau),S)dS=\int_{\partial\Omega^{M-|\sigma|}}e^{i(k+i\tau)\mathcal{L}(S)}b_{0}^{\sigma}((k+i\tau),S)dS,

for a symbol $b_{0}^{\sigma}\in S_{\varrho}^{-|\sigma|}(\partial\Omega^{M-|\sigma|})$ having an asymptotic expansion of the form

\displaystyle b_{0}^{\sigma}(S,k+i\tau)\sim\sum_{p}b_{0,p}^{\sigma}(S)k^{-\sigma-p},

with each $b_{0,p}^{\sigma}$ depending on at most $p$ derivatives of the boundary curvature. From this, an application of the stationary phase lemma gives

\displaystyle\int_{\partial\Omega^{M}}e^{ik\mathcal{L}(S)}a_{0}^{\sigma}((k+i\tau),S)dS\sim\left(\frac{2\pi}{k}\right)^{\frac{M-|\sigma|}{2}}\frac{e^{ik\mathcal{L}(S_{\gamma})}e^{i\pi\text{sgn \,}\partial^{2}\mathcal{L}(S_{\gamma})/4}}{\sqrt{|\det\partial^{2}\mathcal{L}(S_{\gamma})|}}\sum_{j}c_{j,\sigma}k^{-j-|\sigma|}

where $c_{j,\sigma}$ has $3j$ inverse Hessians. Having localized $\hat{\rho}$ near $L_{\gamma}$ , the integral is $O(k^{-\infty})$ unless $M-|\sigma|=q$ , in which case the coefficient of any $k^{-j_{0}}$ in the expansion has $3(j_{0}-|\sigma|)$ inverse Hessians. After separating out the prefactor and using the fact that the derivatives of angles at reflection points depend only on the first order curvature jet, we see that the nonzero $\sigma$ terms contribute $3j-3|\sigma|$ inverse Hessians multiplied by a coefficient which depends smoothly on parameters and a priori on the first $6(j-|\sigma|)$ derivatives of the boundary curvature. However, it is shown in [Zel09] that this can be reduced to $2j$ derivatives of boundary curvature, essentially because in 32 and 33, powers of the higher order phase expansion $g^{\mu}$ vanish to order $3\mu$ . It takes $3\mu$ derivatives to remove this zero, which leaves at most $2j-\mu$ derivatives left to act on the phase or amplitude. The amplitude only depends on the $0$ -jet of $\kappa$ and $g$ depends on the $1$ -jet, but only appears when $\mu\neq 0$ . In either case, at most $2j$ derivatives of $\kappa$ emerge.

We now restrict our attention to the case $|\sigma|=0$ . Since $|x(s_{i})-x(s_{i+1})|\leq\text{diam}(\Omega)\leq|\partial\Omega|$ , we have that for $S$ near $S_{\gamma}$ ,

\displaystyle a_{0}(k+i\tau,S)=b_{0}^{0}(S,k+i\tau)=\left(\frac{k+i\tau}{2\pi}\right)^{\frac{q}{2}}\mathcal{L}(S)\prod_{i=1}^{q}\frac{\cos\vartheta_{i}}{|x(s_{i})-x(s_{i+1})|}+\mathcal{R}_{a_{0},0}(k+i\tau,S),

where the remainder satisfies

\displaystyle|\partial_{z}^{\alpha}\partial_{k}^{\beta}\mathcal{R}_{a_{0},0}(k+i\tau,S)|\leq C_{\alpha,\beta,\varrho,|\partial\Omega|,q}\langle k\rangle^{\varrho|\alpha|-|\beta|-\frac{q}{2}-1}

uniformly and with smooth dependence on parameters, including $\varepsilon$ and the full jet of $\kappa$ . Here, we have chosen only leading order terms in the asymptotic expansion of $a_{1}$ . The contributions of $\mathcal{R}_{a_{0},0}$ and $\frac{\partial}{\partial k}A_{\sigma}\in S_{\varrho}^{-1}$ to the $k^{-j}$ term in the stationary phase expansion have submaximal inverse Hessians for essentially the same reason as the $|\sigma|\geq 1$ terms: an asymptotic expansion of either which starts with $k^{-\widetilde{j}}$ ( $\widetilde{j}\geq 1$ ) contributes to the $k^{-j}$ term in stationary phase expansion through the operator $L_{j-\widetilde{j}}$ , which has at most $3(j-\widetilde{j})$ inverse Hessians.

Thus, for the purposes of calculating wave invariants modulo lower order inverse Hessians, we may replace the integral in Proposition 4.10 by

(35)

\displaystyle\left(\frac{k+i\tau}{2\pi}\right)^{\frac{q}{2}}\int_{\partial\Omega^{q}}e^{i(k+i\tau)\mathcal{L}(S)}\mathcal{L}(S)\prod_{i=1}^{q}\frac{\cos\vartheta_{i}}{|x(s_{i})-x(s_{i+1})|}dS,

to which we apply the stationary phase lemma with parameter $\tau\in[0,1]$ . As noted in Section 4, holomorphy of the regularized resolvent trace allows us to take $\tau\to 0^{+}$ . Denote by $\widetilde{a}(\lambda,S)$ the integrand of 35, corresponding to the amplitude with $|\sigma|=0$ and highest powers of $\varepsilon^{-1}$ . To find explicitly the stationary phase coefficients, we apply the Feynman rules. In Corollary 6.6, we saw the terms which appear in each Feynman amplitide. Most labeled diagrams contribute a Feynman amplitude of zero, since the derivatives of $\mathcal{L}$ are only nonzero if $\max\{|i_{1}-i_{2}|,|i_{1}-i_{3}|,|i_{2}-i_{3}|\}\leq 1$ . Hence, the combinatorial constant depends only on the order of the automorphism group of $3$ -regular graphs on $2j$ vertices. Plugging in $\widetilde{a}$ and $\mathcal{L}$ to the formula in Proposition 6.4, and considering a $3$ -regular graph $\mathcal{G}$ , we sum over all labelings to obtain the Feynman amplitude

	$\displaystyle\sum_{\ell}F_{\mathcal{G},\ell}=$	$\displaystyle\sum_{\ell}\widetilde{a}(S_{\gamma})\prod_{\text{edges}}\prod_{\text{vertices}}\left(ik\partial_{i_{1}(\ell),i_{2}(\ell),i_{3}(\ell)}^{3}\mathcal{L}\right)\left(\frac{i}{k}\widetilde{h}_{j_{1}(\ell)j_{2}(\ell)}\right)$
	$\displaystyle=$	$\displaystyle(\pm i)^{j}k^{-j}\left(c_{\gamma}\varepsilon\right)^{-3j}\widetilde{a}(S_{\gamma})\left(\sum_{i_{1},i_{2},i_{3}}h_{i_{1}}h_{i_{2}}h_{i_{3}}\partial_{i_{1},i_{2},i_{3}}^{3}\mathcal{L}(S_{\gamma})\right)^{2j}.$

Putting $M=q$ , $|\sigma|=0$ , dividing by the order of the automorphism groups, summing over all Feynman diagrams and collecting remainder terms finishes the proof. ∎

7. Future directions

As mentioned in the introduction, this paper is part of a series in which we aim to show that the singular support of the wave trace and the length spectrum are inherently different objects. Our subsequent paper will produce cancellations to arbitrarily high orders in the wave trace. It would also be interesting to both derive degenerate asymptotics and create cancellations in:

•

the wave trace for Riemannian manifolds without boundary,
•

the semiclassical trace formula for the density of states of a Schrödinger operator on a Riemannian manifold with or without boundary, and
•

and the semiclassical trace formula for eigenvalues of large graphs.

In the setting of convex planar domains considered in the present paper, the billiard ball map forms a sort of global Poincaré expression for the billiard flow, which consists of broken bicharacteristics. Formally, our result can then be interpreted as an asymptotic expansion of the trace of a certain boundary operator which quantizes the Poincaré map, analogous to the approach in [SZ02] and [ISZ02]. The authors there use an abstract reduction to a Grushin problem and in future work, it would be interesting to compare their method with the layer potential approach we employ here.

8. Acknowledgements

The first and second authors acknowledge the partial support of the ERC Grant #885707. The third author is grateful to IST Austria and the Schrödinger Institute in Vienna for hosting him during much of the writing of this paper. The authors would also like to thank Hamid Hezari for suggesting the Balian-Bloch approach and the late Steven Morris Zelditch for many helpful conversations in the begining of this project. We are grateful to have known Zelditch for many years and are deeply saddened by his passing. His contributions have made a profound impact on many fields in mathematics. His energy, passion and enthusiasm uplifted all those around him. He was a great inspiration to many mathematicians, including ourselves.

References

[ADSK16] Artur Avila, Jacopo De Simoi, and Vadim Kaloshin. An integrable deformation of an ellipse of small eccentricity is an ellipse. Ann. of Math. (2), 184(2):527–558, 2016.
[AM77] K. G. Andersson and R. B. Melrose. The propagation of singularities along gliding rays. Invent. Math., 41(3):197–232, 1977.
[Axe97] Scott Axelrod. Overview and warmup example for perturbation theory with instantons. In Geometry and physics (Aarhus, 1995), volume 184 of Lecture Notes in Pure and Appl. Math., pages 321–338. Dekker, New York, 1997.
[BB71] R. Balian and C. Bloch. Distribution of eigenfrequencies for the wave equation in a finite domain. II. Electromagnetic field. Riemannian spaces. Ann. Physics, 64:271–307, 1971.
[BB72] R. Balian and C. Bloch. Distribution of eigenfrequencies for the wave equation in a finite domain. III. Eigenfrequency density oscillations. Ann. Physics, 69:76–160, 1972.
[BB74] R. Balian and C. Bloch. Errata: “Distribution of eigenfrequencies for the wave equation in a finite domain. I. Three-dimensional problem with smooth boundary surface” (Ann. Physics 60 (1970), 401–447). Ann. Physics, 84:559, 1974.
[Bia93] Misha Bialy. Convex billiards and a theorem by E. Hopf. Math. Z., 214(1):147–154, 1993.
[Cal22] Keagan G. Callis. Absolutely Periodic Billiard Orbits of Arbitrarily High Order. arXiv e-prints, page arXiv:2209.11721, September 2022.
[CdV84] Y. Colin de Verdière. Sur les longueurs des trajectoires périodiques d’un billard. In South Rhone seminar on geometry, III (Lyon, 1983), Travaux en Cours, pages 122–139. Hermann, Paris, 1984.
[Chr17] T. J. Christiansen. A sharp lower bound for a resonance-counting function in even dimensions. Ann. Inst. Fourier (Grenoble), 67(2):579–604, 2017.
[DG75] Johannes J. Duistermaat and Victor W. Guillemin. The spectrum of positive elliptic operators and periodic bicharacteristics. Invent. Math., 29(1):39–79, 1975.
[dSKW17] Jacopo de Simoi, Vadim Kaloshin, and Qiaoling Wei. Dynamical spectral rigidity among $\mathbb{Z}_{2}$ -symmetric strictly convex domains close to a circle (appendix b coauthored with h. hezari). Ann. of Math. (2), 186(1):277–314, 2017. Appendix B coauthored with H. Hezari.
[Dui96] Johannes J. Duistermaat. Fourier integral operators, volume 130 of Progress in Mathematics. Birkhäuser Boston, Inc., Boston, MA, 1996.
[GK80] V. Guillemin and D. Kazhdan. Some inverse spectral results for negatively curved $2$ -manifolds. Topology, 19(3):301–312, 1980.
[GL18] Colin Guillarmou and Thibault Lefeuvre. The marked length spectrum of Anosov manifolds. ArXiv e-prints, June 2018.
[GM79] Victor Guillemin and Richard Melrose. The Poisson summation formula for manifolds with boundary. Adv. in Math., 32(3):204–232, 1979.
[GWW92] Carolyn Gordon, David L. Webb, and Scott Wolpert. One cannot hear the shape of a drum. Bull. Amer. Math. Soc. (N.S.), 27(1):134–138, 1992.
[HKS18] Guan Huang, Vadim Kaloshin, and Alfonso Sorrentino. On the marked length spectrum of generic strictly convex billiard tables. Duke Math. J., 167(1):175–209, 01 2018.
[Hör03] Lars Hörmander. The analysis of linear partial differential operators. I. Classics in Mathematics. Springer-Verlag, Berlin, 2003. Distribution theory and Fourier analysis, Reprint of the second (1990) edition [Springer, Berlin; MR1065993 (91m:35001a)].
[HZ04] Andrew Hassell and Steve Zelditch. Quantum ergodicity of boundary values of eigenfunctions. Comm. Math. Phys., 248(1):119–168, 2004.
[HZ10] Hamid Hezari and Steve Zelditch. Inverse spectral problem for analytic $(\mathbb{Z}/2\mathbb{Z})^{n}$ -symmetric domains in $\mathbb{R}^{n}$ . Geom. Funct. Anal., 20(1):160–191, 2010.
[HZ12] Hamid Hezari and Steve Zelditch. $C^{\infty}$ spectral rigidity of the ellipse. Anal. PDE, 5(5):1105–1132, 2012.
[HZ22] Hamid Hezari and Steve Zelditch. One can hear the shape of ellipses of small eccentricity. Ann. of Math. (2), 196(3):1083–1134, 2022.
[ISZ02] Alexei Iantchenko, Johannes Sjöstrand, and Maciej Zworski. Birkhoff normal forms in semi-classical inverse problems. arXiv preprint math/0201190, 2002.
[Kac66] Mark Kac. Can one hear the shape of a drum? Amer. Math. Monthly, 73(4, part II):1–23, 1966.
[Kat05] Anatole B. Katok. Billiard table as a playground for a mathematician. In Surveys in modern mathematics, volume 321 of London Math. Soc. Lecture Note Ser., pages 216–242. Cambridge Univ. Press, Cambridge, 2005.
[Kov21] Illya Koval. Local strong Birkhoff conjecture and local spectral rigidity of almost every ellipse. arXiv e-prints, page arXiv:2111.12171, November 2021.
[KS18] Vadim Kaloshin and Alfonso Sorrentino. On the local Birkhoff conjecture for convex billiards. Ann. of Math. (2), 188(1):315–380, 2018.
[KT91] Valeriĭ V. Kozlov and Dmitriĭ V. Treshchëv. Billiards, volume 89 of Translations of Mathematical Monographs. American Mathematical Society, Providence, RI, 1991. A genetic introduction to the dynamics of systems with impacts, Translated from the Russian by J. R. Schulenberger.
[KZ18] Vadim Kaloshin and Ke Zhang. Density of convex billiards with rational caustics. Nonlinearity, 31(11):5214–5234, 2018.
[Mil64] J. Milnor. Eigenvalues of the laplace operator on certain manifolds. Proceedings of the National Academy of Sciences, 51(4):542–542, 1964.
[MM82] Shahla Marvizi and Richard Melrose. Spectral invariants of convex planar regions. J. Differential Geom., 17(3):475–502, 1982.
[OPS88a] Brad Osgood, Ralph Phillips, and Peter Sarnak. Compact isospectral sets of plane domains. Proc. Nat. Acad. Sci. U.S.A., 85(15):5359–5361, 1988.
[OPS88b] Brad Osgood, Ralph Phillips, and Peter Sarnak. Compact isospectral sets of surfaces. J. Funct. Anal., 80(1):212–234, 1988.
[OPS88c] Brad Osgood, Ralph Phillips, and Peter Sarnak. Extremals of determinants of Laplacians. J. Funct. Anal., 80(1):148–211, 1988.
[Pop94] Georgi Popov. Invariants of the length spectrum and spectral invariants of planar convex domains. Comm. Math. Phys., 161(2):335–364, 1994.
[PS92] Vesselin M. Petkov and Luchezar N. Stoyanov. Geometry of reflecting rays and inverse spectral problems. Pure and Applied Mathematics (New York). John Wiley & Sons, Ltd., Chichester, 1992.
[PS17] Vesselin M. Petkov and Luchezar N. Stoyanov. Geometry of the generalized geodesic flow and inverse spectral problems. John Wiley & Sons, Ltd., Chichester, second edition, 2017.
[PT11] Georgi Popov and Petar Topalov. On the integral geometry of Liouville billiard tables. Comm. Math. Phys., 303(3):721–759, 2011.
[PT12] Georgi Popov and Peter Topalov. Invariants of isospectral deformations and spectral rigidity. Comm. Partial Differential Equations, 37(3):369–446, 2012.
[Sib04] Karl Friedrich Siburg. The principle of least action in geometry and dynamics. Number 1844. Springer Science & Business Media, 2004.
[Sjö97] J. Sjöstrand. A trace formula and review of some estimates for resonances. In Microlocal analysis and spectral theory (Lucca, 1996), volume 490 of NATO Adv. Sci. Inst. Ser. C: Math. Phys. Sci., pages 377–437. Kluwer Acad. Publ., Dordrecht, 1997.
[SM85] Toshikazu Sunada and Yozō Matsushima. Riemannian coverings and isospectral manifolds. Annals of Mathematics, 121:169, 1985.
[SZ02] Johannes Sjöstrand and Maciej Zworski. Quantum monodromy and semi-classical trace formulae. Journal de mathématiques pures et appliquées, 81(1):1–33, 2002.
[Tab95] Serge Tabachnikov. Billiards. (No Title), 1995.
[Tay13] Michael Taylor. Partial differential equations II: Qualitative studies of linear equations, volume 116. Springer Science & Business Media, 2013.
[Vig80] Marie-France Vignéras. Variétés riemanniennes isospectrales et non isométriques. Ann. of Math. (2), 112(1):21–32, 1980.
[Vig20] Amir B. Vig. The Inverse Spectral Problem for Convex Planar Domains. ProQuest LLC, Ann Arbor, MI, 2020. Thesis (Ph.D.)–University of California, Irvine.
[Vig21] Amir Vig. Robin spectral rigidity of the ellipse. J. Geom. Anal., 31(3):2238–2295, 2021.
[Vig23a] Amir Vig. Compactness of marked length isospectral sets of birkhoff billiard tables. arXiv preprint arXiv:2310.05426, 2023.
[Vig23b] Amir Vig. The wave trace and birkhoff billiards. Journal of Spectral Theory, 12(3):877–938, 2023.
[WYZ23] Jared Wunsch, Mengxuan Yang, and Yuzhou Zou. The Morse index theorem for mechanical systems with reflections. arXiv e-prints, page arXiv:2308.16162, August 2023.
[Zel98] Steve Zelditch. The inverse spectral problem for surfaces of revolution. J. Differential Geom., 49(2):207–264, 1998.
[Zel00] S. Zelditch. Spectral determination of analytic bi-axisymmetric plane domains. Geom. Funct. Anal., 10(3):628–677, 2000.
[Zel04a] Steve Zelditch. Inverse resonance problem for $\mathbb{Z}_{2}$ -symmetric analytic obstacles in the plane. In Christopher B. Croke, Michael S. Vogelius, Gunther Uhlmann, and Irena Lasiecka, editors, Geometric Methods in Inverse Problems and PDE Control, pages 289–321, New York, NY, 2004. Springer New York.
[Zel04b] Steve Zelditch. The inverse spectral problem. In Surveys in differential geometry. Vol. IX, volume 9 of Surv. Differ. Geom., pages 401–467. Int. Press, Somerville, MA, 2004. With an appendix by Johannes Sjöstrand and Maciej Zworski.
[Zel04c] Steve Zelditch. Inverse spectral problem for analytic domains. I. Balian-Bloch trace formula. Comm. Math. Phys., 248(2):357–407, 2004.
[Zel09] Steve Zelditch. Inverse spectral problem for analytic domains. II. $\mathbb{Z}_{2}$ -symmetric domains. Ann. of Math. (2), 170(1):205–269, 2009.
[Zwo98] Maciej Zworski. Poisson formula for resonances in even dimensions. Asian J. Math., 2(3):609–617, 1998.