Eigenfunction asymptotics and spectral Rigidity of the ellipse

This paper is part of a series concerning the isospectral problem for an ellipse. In this paper, we study Cauchy data of eigenfunctions of the ellipse with Dirichlet or Neumann boundary conditions. Using many classical results on ellipse eigenfunctions, we determine the microlocal defect measures of the Cauchy data of the eigenfunctions. The ellipse has integrable billiards, i.e. the boundary phase space is foliated by invariant curves of the billiard map. We prove that, for any invariant curve $C$, there exists a sequence of eigenfunctions whose Cauchy data concentrates on $C$. We use this result to give a new proof that ellipses are infinitesimally spectrally rigid among $C^{\infty}$ domains with the symmetries of the ellipse.

This note is part of a series [HeZe12,HeZe19] on the the inverse spectral problem for elliptical domains E ⊂ R 2 . In [HeZe12], it is shown, roughly speaking, that an isospectral deformation of an ellipse through smooth domains (but not necessarily real analytic) which preserves the Z 2 × Z 2 symmetry is trivial. In [HeZe19] it is shown that ellipses of small eccentricity are uniquely determined by their Dirichlet (or, Neumann) spectra among all C ∞ domains, with no analyticity or symmetry assumptions imposed. In both [HeZe12,HeZe19], the main spectral tool is the wave trace singularity expansion and the special form it takes in the case of ellipses. In this article, we take the dual approach of studying the asymptotic concentration in the phase space B * ∂E of the Cauchy data u b j of Dirichlet (or, Neumann) eigenfunctions u j of elliptical domains in the unit coball bundle of the boundary ∂E. In Theorem 1, we show that, for every regular rotation number of the billiard map in the 'twist interval', there exists a sequence of eigenfunctions whose Cauchy data concentrates on the invariant curve with that rotation number in B * ∂E. The proof uses the classical separation of variables and one dimensional WKB analysis.
We are interested here in the quantum limits of the Cauchy data (1) of an orthonormal basis of eigenfunctions of an ellipse, i.e. in the asymptotic limits of the matrix elements of zeroth order semi-classical pseudo-differential operators Op (a) on ∂E with respect to the L 2normalized Cauchy data of eigenfunctions. We note that ρ b j is normalized so that ρ b j (I) = 1 and is a positive linear functional, hence all possible weak* limits are probability measures on the unit coball bundle B * ∂Ω. Moreover, ρ j (N (λ) * Op (a)N (λ)) = ρ j (Op (a)), so that the quantum limits are quasi-invariant under the billiard map (see [HaZe04] for precise statements). In Theorem 1 we determine the quantum limits of sequences in (2) for an ellipse. The proof uses many of the prior results on WKB formulae for ellipse eigenfunctions, especially those of [KeRu60,WaWiDu97,Sie97].
In large part, our interest in matrix elements (2) owes to the fact that the Hadamard variational formulae for eigenvalues of the Laplacian with Dirichlet boundary condition expresses the eigenvalue variations as the special matrix elements (2) given by, of the domain variationρ (not to be confused with ρ j ) against squares of the Cauchy data (see Section 5.1). As stated in Corollary 2, the limits of such integrals over all possible subsequences of eigenfunctions determines the 'Radon transform' ofρ over all possible invariant curves for the billiard map. Under an infinitesimal isospectral deformation, all of the limits are zero. We use this result to give a new proof of the spectral rigidity result in [HeZe12]; see Theorem 4 and Corollary 5. The principal motivation for studying the inverse Laplace spectral problems for ellipses stems from the Birkhoff conjecture that ellipses are the only bounded plane domains with completely integrable billiards. Strong recent results, due to A. Avila, J. de Simoi, V. Kaloshin, and A. Sorrentino [AvdSKa16,KaSo18] have proved local versions of the Birkhoff conjecture using a weaker notion of integrability known as 'rational integrability', i.e. that periodic orbits come in oneparameter families, namely invariant curves of the billiard map with rational rotation number. In this article, Bohr-Sommerfeld invariant curves play the principal role rather than curves of periodic orbits.
1.1. Statement of results. The first result pertains to concentration of Cauchy data of sequences ϕ j of Dirichlet (resp. Neumann) eigenfunctions on invariant curves of the billiard map of an ellipse. We denote E by x 2 a 2 + y 2 b 2 ≤ 1, 0 ≤ b < a, and choose the elliptical coordinates (ρ, ϑ) by (x, y) = (c cosh ρ cos ϑ, c sinh ρ sin ϑ). Here, We denote the angular Hamiltonian, which we will also call the action, by I = p 2 ϑ /c 2 + cos 2 ϑ. The invariant curves of β are the level sets of I. The range of I is called the action interval. There is a natural measure dµ α on each level set I = α called the Leray measure which is invariant under β and the flow of I. We refer to Section 2 for detailed definitions and properties involving the billiard map of an ellipse, actions, invariant curves, and the Leray measure.
Theorem 1. Let E be an ellipse. For any α in the action interval of the billiard map of E, there exists a sequence of separable (in elliptical coordinates) eigenfunctions {ϕ j } of eigenvalue λ 2 j whose Cauchy data concentrates on the level set {I = α}, in the sense that, for any zeroth order semi-classical pseudo-differential operator Op (a) on B * ∂E with principal symbol a 0 , In particular, Corollary 2. In the special case when the symbol a(ϑ, p ϑ ) =ρ(ϑ) is only a function of the base variable ϑ, where ds = c 2 (cosh 2 ρ max − cos 2 ϑ) dϑ is the arclength measure.
Remark 3. If we denote η to be the symplectic dual variable of the arclength s, then our quantum limit can be expressed as For the proof, see our computation of 1 − |η| 2 in the proof of Corollary 18. The appearance of the (non-invariant) factors 1 − |η| 2 and 1/ 1 − |η| 2 is consistent with the result of [HaZe04], where the quantum limits of boundary traces of ergodic billiard tables are studied.
To our knowledge, Theorem 1 is the first result on microlocal defect measures of Cauchy data of eigenfunctions in non-ergodic cases. See Section 1.3 for related results. One of the difficulties in determining the limits of (2) is that the Cauchy data u b j are not L 2 normalized. It is shown in [HaTa02, Theorem 1.1] that there exists C, c > 0 so that c ≤ ||λ −1 j ∂ ν ϕ j || L 2 (∂Ω) ≤ C for Dirichlet eigenfunctions of Euclidean plane domains (and more general non-trapping cases). Hence the L 2 normalization in (2) is rather mild. On the other hand, the corresponding inequalities do not hold in general for Neumann eigenfunctions. As pointed out in [HaTa02,Example 7], there are simple counter-examples to any constant upper bound on the unit disc (whispering gallery modes). There do exist positive lower bounds for convex Euclidean domains. Hence, in the case of an ellipse, the L 2 normalization in (2) is necessary to obtain limits.
1.2. Spectral rigidity. Before stating the results, we review the main definitions. An isospectral deformation of a plane domain Ω 0 is a one-parameter family Ω t of plane domains for which the spectrum of the Euclidean Dirichlet (or Neumann, or Robin) Laplacian ∆ t is constant (including multiplicities). The deformation is said to be a C 1 deformation through C ∞ domains if each Ω t is a C ∞ domain and the map t → Ω t is C 1 . We parameterize the boundary ∂Ω t as the image under the map where ρ t ∈ C 1 ([0, t 0 ], C ∞ (∂Ω)). The first variation is defined to beρ(x) := d dt | t=0 ρ t (x). An isospectral deformation is said to be trivial if Ω t = Ω 0 (up to isometry) for sufficiently small t. A domain Ω 0 is said to be spectrally rigid if all C ∞ isospectral deformations are trivial.
In [HeZe12] the authors proved a somewhat weaker form of spectral rigidity for ellipses, with 'flatness' replacing 'triviality'. Its main result is the infinitesimal spectral rigidity of ellipses among C ∞ plane domains with the symmetries of an ellipse. We orient the domains so that the symmetry axes are the x-y axes. The symmetry assumption is then that ρ t is invariant under (x, y) → (±x, ±y). The variation is called infinitesimally spectrally rigid ifρ 0 = 0.
The main result of [HeZe12] is: Theorem 4. Suppose that Ω 0 is an ellipse, and that Ω t is a C 1 Dirichlet (or Neumann) isospectral deformation of Ω 0 through C ∞ domains with Z 2 × Z 2 symmetry. Let ρ t be as in (6). Thenρ = 0.
Corollary 5. Suppose that Ω 0 is an ellipse, and that t → Ω t is a C 1 Dirichlet (or Neumann) isospectral deformation through Z 2 × Z 2 symmetric C ∞ domains. Then ρ t must be flat at t = 0.
The proof of Theorem 4 in [HeZe12] used the variation of the wave trace. In the original posting (arXiv:1007.1741) the authors used a more classical Hadamard variational formula for variations of individual eigenvalues λ j (t), which appears in Section 5.1. The authors rejected this approach in favor of the one appearing in [HeZe12] because it was thought that this argument was invalid when the eigenvalues were multiple. When a multiple eigenvalue of a 1-parameter family L t of operators is perturbed, it splits into a collection of branches which in general are not differentiable in t. Moreover, the authors assumed that the variational formula would express the variation in terms of special separable eigenfunctions (see Section 3). This created doubt that one could use the variational formula for individual eigenvalues. Instead, the authors used the variational formula for the trace of the wave group or equivalently for spectral projections, which are symmetric sums over all of the branches into which an eigenvalue splits.
However, as we show in this article, the original variational formulae were in fact correct even in the presence of multiplicities. The first point is that the non-differentiability issue does not arise for an isospectral deformation since no splitting occurs. Second, the vanishing of the variation of eigenvalues implies that the infinitesimal variationρ is orthogonal to squares of all (Dirichlet) eigenfunctions in the eigenspace, and in particular the separable ones. More precisely, we prove that ∂Eρ |u b j | 2 ds = 0.
Then by Corollary 5, we obtain that for every α in the action interval one has (7) I=αρ dν α = 0.
In the final step we calculate the measure dν α and provide two proofs, one via inverting an Abel transform and another using the Stone-Weierstrass theorem, that (7) impliesρ = 0. The proof in the Neumann case is similar and will be provided.

1.3.
Related results and open problems. Quantum limits of Cauchy data on manifolds with boundary have been studied in [HaZe04,ChToZe13] in the case where the billiard map β is ergodic. To our knowledge, they have not been studied before in non-ergodic cases. Theorem 1 shows that, as expected, Cauchy data of eigenfunctions localize on invariant curves for the billiard map rather than delocalize as in ergodic cases. L 2 norms of Cauchy data of eigenfunctions are studied in [HaTa02] in the Dirichlet case and in [BaHaTa18] in the Neumann case. Further results on the quasi-orthonormality properties of Cauchy are studied in [BFSS02,HHHZ15].
The study of eigenfunctions in ellipses has a long literature and we make substantial use of it. In particular, we quote several articles in the physics literature, in particular [WaWiDu97,Sie97], and several in mathematics [KeRu60,BaBu91], for detailed analyses of eigenfunctions of the quantum ellipse. There is also a series of articles of G. Popov and P. Topalov (see e.g. [PoTo03,PT16]) on the use of KAM quasi-modes to study Laplace inverse spectral problems. In particular, in [PT16], Popov-Topalov also give a new proof of the rigidity result of [HeZe12] and extend it to other settings. The approach in this article is closely related to theirs, although it does not seem that the authors directly studied Cauchy data of eigenfunctions of an ellipse.
The multiplicity of Laplace eigenvalues of an ellipse appears to be largely an open problem. It is a non-trivial result of C.L. Siegel that the multiplicities are either 1 or 2 in the case of circular billiards; multiplicity 1 occurs for, and only for, rotationally invariant eigenfunctions. The Laplacians of the family of ellipses x 2 a 2 + y 2 b 2 = 1 form an analytic family containing the disk Laplacian, and one might try to use analytic perturbation theory to prove the following, Conjecture 6. For a generic class of ellipses the multiplicity of each eigenvalue is ≤ 2.
1.4. Quantum Birkhoff conjecture. As mentioned above, ellipses have completely integrable billiards, and the classical Birkhoff conjecture is that elliptical billiards are the only completely integrable Euclidean billiards with convex bounded smooth domains. Despite much recent progress, the Birkhoff conjecture remains open.
The eigenvalue problem on a Euclidean domain is often called 'quantum billiards' in the physics literature (see e.g. [WaWiDu97]). One could formulate quantum analogues of the Birkhoff conjecture in several related but different ways. The quantum analogue of the Birkhoff conjecture is presumably that ellipses are the only 'quantum integrable' billiard tables. A standard notion of quantum integrability is that the Laplacian commutes with a second, independent, (pseudodifferential) operator; we refer to [ToZe03] for background on quantum integrability. In Section 3, we explain that the ellipse is quantum integrable in that one may construct two commuting Schrödinger operators with the same eigenfunctions and eigenvalues. The symbol of the second operator then Poisson commutes with the symbol of the Laplacian, hence the billiard dynamics and billiard map are integrable. A related version is that one can separate variables in solving the Laplace eigenvalue problem. It is not obvious that these two notions are equivalent; in Section 3 we use both separation of variables and existence of commuting operators in studying the ellipse. Classical studies of separation of variables and its relation to integrability go back to C. Jacobi, P. Stäckel, L. Eisenhart and others, and E.K. Sklyanin has studied the problem more recently. We do not make use of their results here.
Quantum integrability is much stronger than classical integrability, and one might guess that it is simpler to prove the quantum Birkhoff conjecture than the classical one. Wave trace techniques as in [HeZe12,HeZe19] reduce Laplace spectral determination and rigidity problems to dynamical inverse or rigidity results. The wave trace only 'sees' periodic orbits and is therefore well-adapted to results on rational integrability. The dual approach through eigenfunctions studied in this article gives a different path to the quantum Birkhoff conjecture, in which rational integrability and periodic orbits play no role.
Acknowledgment. We thank Luc Hillairet for a discussion which prompted the revival of this note.

Classical billiard dynamics
In this, and the next, section, we review some background definitions and results on the classical and quantum elliptical billiard. We follow the notation of [Sie97]; see also [BaBu91,WaWiDu97].
An ellipse E is a plane domain defined by, Here, a, resp. b, is the length of the semi-major (resp. semi-minor) axis. The ellipse has foci at (±c, 0) with c = √ a 2 − b 2 and its eccentricity is e = c a . Its area is πab, which is fixed under an 6 HAMID HEZARI AND STEVE ZELDITCH isospectral deformation. We define elliptical coordinates (ρ, ϑ) by (x, y) = (c cosh ρ cos ϑ, c sinh ρ sin ϑ).
The coordinates are orthogonal. The lines ρ = constant are confocal ellipses and the lines ϑ = constant are confocal hyperbolas. In the special case of the disc, we have c = 0, but we assume henceforth that c = 0.
2.1. Action variables for the billiard flow. The billiard flow on the ellipse E is the (broken) geodesic flow of the Hamiltonian H = p 2 x + p 2 y on T * E, which follows straight lines inside E and reflects on ∂E according to equal angle law of reflection.
Action-angle variables on T * E are symplectic coordinates in which the billiard flow of the ellipse is given by Kronecker flows on the invariant Lagrangian submanifolds. We refer to [Ar89] for the general principles and to [Sie97] for the special case of the ellipse. Let p ρ and p ϑ be the symplectic dual variables corresponding to the elliptic coordinates ρ and ϑ, respectively. The two conserved quantities of the system are the energy (the Hamiltonian) H and the angular Hamiltonian I (which we also call the action), given in the coordinates (ρ, p ρ , ϑ, p ϑ ), by In the notation of [Tab97], I = cos 2 θ cosh 2 ρ + sin 2 θ cos 2 ϑ, where θ is the angle between a trajectory of the billiard flow and a tangent vector to the confocal ellipse with parameter ρ. Note also that by the notation of [Sie97], I = 1 + L 1 L 2 c 2 H where L 1 L 2 is the product of two angular momenta about the two foci. The values of I are restricted to 0 ≤ I ≤ a 2 c 2 = cosh 2 (ρ max ). The upper limit I = cosh 2 (ρ max ) corresponds to the motion along the boundary and the lower limit I = 0 corresponds to the motion along the minor axis. Moreover, there are two different kinds of motion in the ellipse depending on the sign of I. For 1 < I < cosh 2 (ρ max ) the trajectories have a caustic in the form of a confocal ellipse. For 0 < I < 1 the caustic of the motion is a confocal hyperbola and the trajectories cross the x-axis between the two focal points. Both kinds of motions are separated by a separatrix which consists of orbits with I = 1 that go through the focal points of the ellipse.
In terms of H and I, the canonical momenta, are given by

Therefore, the action variables are
In fact these are the actions for the half-ellipse 0 ≤ ϕ ≤ π. The integrals can be calculated in terms of I using elliptic integrals of first and second kind (See [Sie97]). The actions will play a key role in Section 3.4 in the description of Bohr-Sommerfeld quantization conditions for the eigenvalues of the Laplacian. that go through the focal points of the ellipse. Figure 1 shows a Poincaré section through the classical motion in an ellipse for different initial conditions. The boundary of the billiard is chosen as the surface of the section and on this surface the reflections are described in Birkhoff coordinates: s is the arclength along the billiard boundary and p is the cosine of the angle between the outgoing trajectory and the tangent to the boundary at the reflection point. The two lines through (s, p) = (0, 0) mark the separatrix of the motion. The lines inside the separatrix correspond to the librational motion with ↵ < 0 and the lines outside the separatrix correspond to the rotational motion  2.2. Billiard map, invariant curves, Leray measure, and action-angle variables. The billiard map of an ellipse E (or in general any smooth domain) is a cross section to the the billiard flow on S * ∂E E, which we always identify with B * ∂E and call it the phase space of the boundary. To be precise, the billiard map β is defined on B * ∂E as follows: given (s, η) ∈ T * ∂E, with s being the arc-length variable measured in the counter-clockwise direction from a fixed point say s 0 , and |η| ≤ 1, we let (s, ζ) ∈ S * E be the unique inward-pointing unit covector at s which projects to (s, η) under the map T * ∂E E → T * ∂E. Then we follow the geodesic (straight line) determined by (s, ζ) to the first place it intersects the boundary again; let s ∈ ∂E denote this first intersection. (If |η| = 1, then we let s = s.) Denoting the inward unit normal vector at s by ν s , we let ζ = ζ + 2(ζ · ν s )ν s be the direction of the geodesic after elastic reflection at s , and let η be the projection of ζ to T * s Y . Then we define β(s, η) = (s , η ). A theorem of Birkhoff asserts that billiard map preserves the natural symplectic form ds ∧ dη on B * ∂E, i.e. β * (ds ∧ dη) = ds ∧ dη. In the literature, the coordinates (s, θ) are commonly used for phase space of the boundary, where θ ∈ [0, π] is the angle that ζ makes with the positive tangent direction at s. In these coordinates, An invariant curve is a curve (connected or not) on the phase space that is invariant. The phase space B * ∂E of the ellipse E is in fact foliated with invariant curves. More precisely, Although I ϑ is the classical angular action on B * ∂E, but we shall call I the action as it is more convenient and is related to I ϑ via the one-to-one correspondence (10). As is evident from the Figure 1, the separatrix curve I = 1 divides the phase space into two types of open sets, the exterior corresponding to trajectories with confocal elliptical caustics (1 < I < cosh 2 ρ max ) and the interior to trajectories with confocal hyperbolic caustics (0 < I < 1).
2.2.1. Leray measure. On each level set I = α of I, there is a natural measure dµ α called the Leray measure which in invariant under β and the flow generated by I. In the symplectic coordinates (ϑ, p ϑ ), and on I = α, it is given by Since dϑ ∧ dI = ∂I ∂p ϑ dϑ ∧ dp ϑ , we obtain that Here, x + = x if x > 0 and is zero otherwise. Up to a scalar multiplication, dµ α is a unique measure that is invariant under β and the flow of I.

2.2.2.
Action-angle variables and rotation number. The billiard map has a Birkhoff normal form around each invariant curve in B * ∂E. That is, in the symplectically dual angle variable ι to I, the billiard map has the form, β(I, ι) = (I, ι + r(I)), where r is often called the rotation number of the invariant curve. An explicit formula is given for it in [Tab97] (3.5), [CaRa10] (section 4.3 (11)) and [Ko85]. Then, if 0 < I < 1, Also, if 1 < I < cosh 2 (ρ max ) then Definition: We define the range of the action variable I as the action interval, i.e. the interval [0, cosh 2 (ρ max )], and the range of r(I) as the rotation interval.

Quantum elliptical billiard
The Helmholtz equation in elliptical coordinates takes the form, The quantum integrability of ∆ owes to the fact that this equation is separable. We put where = λ −1 and α is the separation constant. Here, 'PBC' stands for 'periodic boundary conditions'; DBC (resp. NBC) stands for Dirichlet (resp. Neumann) boundary conditions. Thus, we consider pairs ( , α) where there exists a smooth solution of the two boundary problems.
Each of the angular and radial equations above is an eigenvalue problem for a semiclassical Schrödinger operator with boundary conditions on a finite interval. These commuting operators are given by Op (I); I = p 2 ϑ /c 2 + cos 2 (ϑ).
As G(−ϑ) is a solution whenever G(ϑ) is, we restrict our attention to 2π-periodic solutions to the angular equation which are either even or odd. One can then see that: Remark 8. In order to obtain solutions well-defined on the line segment joining the foci, i.e. at ρ = 0, solutions to the radial equation must satisfy the boundary condition F (0) = 0 in case the solution G is even and F (0) = 0 in case G is odd. In these cases the solutions F are also respectively even and odd functions.

Mathieu and modified Mathieu characteristic numbers.
For each fixed , the angular problem is a Sturm-Liouville problem and thus there exist real valued sequences {a n ( )} ∞ n=0 and {b n ( )} ∞ n=1 so that it has 2π-periodic non-trivial solutions -even solutions if α = a n ( ) and odd solutions if α = b n ( ). Here even or odd is with respect to ϑ → −ϑ, or equivalently y → −y. We represent the corresponding solutions by G e n (ϑ, ) and G o n (ϑ, ), respectively. The even indices correspond to π-periodic solutions, thus they must be invariant under ϑ → π − ϑ, or equivalently be even with respect to x → −x. Solutions with odd indices have anti-period π and correspond to odd solutions in the x variable. The sequences a n ( ) and b n ( ) are related to the standard Mathieu characteristic numbers of integer orders a n (q) and b n (q) by (18) a n ( ) = 1 2 + a n (q) 4q , b n ( ) = 1 2 + b n (q) 4q , q = c 2 4 2 .
Thus using the wellknown properties of a n and b n , for > 0 we have The existence of the point of intersection of the curves a n ( ) with A m ( ), and b n ( ) with B m ( ) are guaranteed by: Theorem 9 (Neves [Ne10]). For each (m, n), there is a unique positive solution q to each of the equations a n (q) = A m (q) and b n (q) = B m (q).
Hence the same statement holds for the equations (22) by the correspondence (18). The frequencies λ j of E are obtained by sorting {λ e mn , λ o mn ; (m, n) ∈ N 2 } in increasing order. 3.3. Symmetries classes. The irreducible representations of the Z 2 × Z 2 symmetry group are real one-dimensional spaces, so that there exists an orthonormal basis of eigenfunctions of the ellipse which are even or odd with respect to each Z 2 symmetry, i.e. have one of the four symmetries (even, even), (even, odd), (odd, even), (odd, odd), where the first and the second entries correspond to symmetries with respect to x → −x and y → −y, respectively. Given the above discussion the symmetric eigenfunctions are:  Figure 2 shows the symmetries classes of eigenfunctions distinguished by their probability densities. It is possible that two symmetric eigenfunctions correspond to the same eigenvalue, or it is possible that they correspond to different eigenvalues.

Semiclassical actions and Bohr-Sommerfeld quantization conditions for the ellipse.
Graphs of the one-dimensional classical potentials are given in [WaWiDu97, Figure 1]. The potential − cosh 2 ρ for Op (J) in (15) is a potential barrier with a single local maximim which is symmetric around the vertical line through the local maximum. The classical potential cos 2 ϑ underlying Op (I) in (16) is a double-well potential on the circle. Thus, there exists a separatrix curve corresponding to the two local maxima of the potential, which divides the two-dimensional phase space into two regions. Inside the phase space curve, the level sets of the potential are 'circles' paired by the left right symmetry across the vertical line through the local maximum at π. Outside the separatrix, the level sets have non-singular projections to the base, i.e. are roughly horizontal. It is more important for our purposes to determine the lattice of semi-classical eigenvalues in terms of classical and quantum action variables. The WKB (or EKB) quantization for the actions are given in [Sie97,(33)] (see also [KeRu60] for the original reference). Up to O( 2 ) terms they have the form: , I ϑ = (n + 1 2 ) , m, n = 0, 1, 2, . . . . There is a discontinuity at I = 1 due to the separatrix curve, but it is not important for our problem and we ignore it.
3.4.1. Semiclassical action. In fact for each of the eight Bohr-Sommerfeld quantization condition above there is a version which is valid to all orders in which are essentially given by the quantum Birkhoff normal form around each orbit under consideration. To be precise, there exist eight so called semiclassical actions S e/o,1 ± ,ρ/ϑ (α), where the choices of e/o corresponds to even or odd in the y (equivalently in the ϑ) variable, of 1 + or 1 − to I > 1 or I < 1, and ρ or ϑ to actions in the ρ or ϑ variable, respectively. Each of the eight semiclassical actions has an asymptotic expansion of the form where S 0 (α) is the corresponding classical action which is I ρ | I=α for S e/o,1 ± ,ρ and I ϑ | I=α for where > 0 is arbitrary, however the remainder estimates in the asymptotic expansions depend on . There are versions of BSQC in the literature that are valid uniformly near the separatrix but we do not need it here. We also point out that the Maslov indices are not ignored but absorbed in the corresponding subleading terms S 1 (α). 3.5. Keller-Rubinow algorithm. In this section we explore the procedure of finding e mn corresponding to eigenvalues associated to invariant curves outside the separatrix (i.e. 1 + case) whose eigenfunctions are even in the ϑ variable. All other cases follow a similar procedure and we shall drop the superscripts for convenience.
We are in search of solutions to equation (28) which, in our convenient notation, are given by S ρ (α ρ m ( )) = m , S ϑ (α ϑ n ( )) = n , respectively. Following [KeRu60], we divide these two equations to obtain, The expression A (α) has a classical expansion with principal term which is a positive monotonic function on the interval [1, cosh 2 ρ max ] (See [KeRu60], page 41). Hence, if we choose r in the range of A 0 (α) on the domain [1 + 2 , cosh 2 ρ max − 2 ], then for sufficiently small there is a unique solution α to the equation A (α) = r in the slightly larger interval [1 + , cosh 2 ρ max − ], accepting an expansion of the form: It is manifestly the the inverse function of A h (α) and its formal power series coefficients α (k) (r) are smooth functions of r. The principal term α (0) is the inverse function of A 0 (α). By this definition, the solution to (31) is α( , m/n) whenever m/n belongs to A 0 [1 + 2 , cosh 2 ρ max − 2 ], which is a bounded closed interval in (0, ∞). In particular m/n is bounded above and below by positive constants K 1 and K 2 : This is the eligible sector of lattice points for our eigenvalue problem outside the separatrix. Plugging α( , m/n) into the angular BSQC, i.e. the second equation of (30), (the radial one follows immediately from the angular one and (31)), we arrive at the quantization condition for the eigenvalues of E: (35) Q( , m, n) := 1 n S ϑ (α( , m/n)) = .
We claim that for m and n sufficiently large, this equation has a unique solution mn in a sufficiently small interval [0, 0 ], or equivalently the function Q(·, m, n) has a unique fixed point. Now, since for 0 sufficiently small, and n sufficiently large Q(·, m, n) maps [0, 0 ] into itself and ∂Q ∂ ( , m, n) < 1 2 in this interval. The claim follows by the Banach contraction principle.
Remark 11. Since there are many functions α used, it is important to highlight their relations and differences. If we evaluate α ( , r) defined in (33), at = m,n and r = m n , we get the common value of (29). In short, α mn , m n = α ρ m ( mn ) = α ϑ n ( mn ). We also note that the function α (0) (r), with parentheses around 0, is the principal term of α( , r) and should not be confused with α ρ 0 ( ) or α ϑ 0 ( ). In fact, the above procedure provides an asymptotic for λ mn = 1/ mn and gives a sharper result than previously known: Proposition 12. The frequencies λ e/o mn of E associated to invariant curves outside the separatrix curve, and away from it, correspond to lattice points (m, n) ∈ N 2 in the sector and satisfy the asymptotic property, The same asymptotic formula holds for the frequencies λ e/o mn associated to invariant curves inside the separatrix curve, except in this case the sector of lattice points is: The effects of even/odd are only reflected in the remainder term O(1/n), which in addition depends on the distance from the separatix. Note that the explicit formulas for I ϑ and I ρ (hence for α (0) ) in terms of elliptic integrals are different for the inside and outside the separatrix curve (See for example [Sie97]).  (23) along 'ladders' or 'rays' in the action lattice (m, n) ∈ N 2 . In particular, different rays correspond to different invariant Lagrangrian submanifolds for the billiard flow. It is simpler to use the billiard map and then to relate rays in the joint spectrum to invariant curves for the billiard map. Given an invariant curve, inside or outside the separatrix, we wish to find a ray in the joint spectrum for which the associated eigenfunctions concentrate on the curve. Since the WKB method is highly developed in dimension one, it suffices for our purposes to locate the ray in N 2 which corresponds to the invariant curve. The corresponding eigenfunctions will then concentrate on the corresponding Lagrangian submanifolds. We shall only prove (1), as the proof of (2) is similar. Furthermore, we shall only focus on the even case because the proof for the odd case is identical. Fix α ∈ (1, cosh 2 ρ max ). We choose > 0 so that α ∈ [1 + 2 , cosh 2 ρ max − 2 ]. Let mn be the sequence we found in Section 3.5 associated to the level curves outside the separatrix and to even eigenfunctions (even in the y variable). By Remark 11, it suffices to show that there is a subsequence (m j , n j ) along which We choose r 0 by α (0) (r 0 ) = α (recall that α (0) is monotonic) and choose a sequence of lattice points (m j , n j ) ∈ N 2 in the eligible sector (34) such that , which proves Theorem 1 in the Dirichlet case. The Neumann case is essentially the same; we omit the details.

Hadamard variational formulae for isospectral deformations
We consider the Dirichlet (resp. Neumann) eigenvalue problems for a one parameter family of Euclidean plane domain Ω t , where Ω 0 = E is an ellipse: Here, ∂ νt is the interior unit normal to ∂Ω t . When λ 2 j (0) is a simple eigenvalue, then under a C 1 deformation the eigenvalue moves in a C 1 curve λ 2 j (t). When λ 2 j (0) is a multiple eigenvalue, then in general the eigenvalue may split into branches. Examples in [Ka95] show that eigenfunctions do not necessarily deform nicely if the deformation is not analytic. Hence we cannot even assume that eigenfunctions are C 1 if the deformation is only C 1 . However, we assume in this section that the deformation is isospectral. In this case, a multiple eigenvalue does not change multiplicity under the deformation, and therefore there is no splitting into branches.
When an eigenvalue has multiplicity > 1, there exists an orthonormal basis (known as the Kato-Rellich basis) of the eigenspace which moves smoothly under the deformation. The multiple eigenvalue splits under a generic perturbation and one can only expect a perturbation formula along each path. When we assume that the deformation is isospectral, hence that the eigenvalue does not split (or even change) along the deformation, then there exists a Kato-Rellich basis for the eigenspace. 5.1. Hadamard variational formulae. As in the introduction, we parameterize the deformation by a function ρ t on ∂E so that ∂Ω t is the graph of ρ t over ∂Ω 0 = ∂E in the sense that ∂Ω t = {x + ρ t (x)ν x : x ∈ ∂Ω 0 }. Ifρ := d dt ρ t | t=0 = 0, then the first order variation of eigenvalues is the same as for the deformation by x + tρ(x)ν x . In this section we review the Hadamard variational formula in the case of simple eigenvalues. We refer to [HeZe12, Section 1] for background on the Hadamard variational formula.
When λ 2 j (0) is a simple eigenvalue (i.e. of multiplicity one) with L 2 -normalized eigenfunction ϕ j , then Hadamard's variational formula for plane domains is that where ds is the induced arc-length measure. Hence, under an infinitesimal isospectral deformation we have, for every simple eigenvalue, Dirichlet: Hadamard's variational formula is actually a variational formula for the variation of the Green's functions G(λ, x, y) with the given boundary conditions. In the Dirichlet case it states thaṫ G(λ, x, y) = − ∂Ω 0 ∂ ∂ν 1 G(λ, q, x) ∂ ∂ν 1 G(λ, q, y)ρds.
The formula (39) follows if we compare the poles of order two on each side. The same comparison shows that if the eigenvalue λ 2 j (0) is repeated with multiplicity m(λ j (0)) and if {λ jk (t)} m(λ j (0)) j=1 is the perturbed set of eigenvalues, then Here is any ONB of the repeated eigenvalue λ 2 j (0). There exist similar Hadamard variational formulae in the Neumann case. When the eigenvalue is simple, we have Neumann: 5.2. Hadamard variational formula for an isospectral deformation. We now assume that the deformation is isospectral. As mentioned above, there exists a Kato-Rellich basis which moves smoothly under the deformation. In fact, we show that for an isospectral deformation every eigenfunction has a smooth deformation along the path. In the following −∆ t denotes the Dirichlet (resp. Neumann) Laplacian on Ω t .
Proof. Let λ 2 j (0) be the eigenvalue of ϕ j (0), of multiplicity m j ≥ 1, and γ be a circle in C centered at λ 2 j (0) such that no other eigenvalues of −∆ 0 are in the interior of γ or on γ. We define where R t (z) = (−∆ t − z) −1 is the resolvent of −∆ t . By the Cauchy integral formula, it is clear that P 0 is the orthogonal projector onto the eigenspace of λ 2 j (0). Since the eigenvalues {λ 2 j,k (t)} m j k=1 vary continuously in t, for t small these are the only eigenvalues of −∆ t in γ. Therefore, in general, P t is the total projector (the direct sum of projectors) associated with {λ 2 j,k (t)} m k=1 . The operator P t is C 1 in t, since the resolvent (hence, Green's function) is C 1 in t (see [Ka95,Theorem II.5.4]). Now assume Ω t is an isospectral deformation. Since the spectrum is constant along the deformation, P t projects every function on Ω t onto an eigenfunction of Ω t of eigenvalue λ 2 j (0). Let f t be a C 1 family of smooth diffeomorphisms from Ω t to Ω 0 with f 0 = Id. Then must be an eigenfunction of −∆ t of eigenvalue λ 2 j (0).
We are now in position to prove: Lemma 17. Suppose that Ω t is a C 1 isospectral deformation. Then for any eigenfunction ϕ j of Ω 0 , Proof. Let ϕ j (0) be any eigenfunction of Ω 0 and ϕ j (t) be the C 1 deformation of eigenfunction of Ω t provided by Lemma 16. For t > 0, the eigenvalue problem for the isospectral deformation is pulled back to Ω 0 by a C 1 family diffeomorphisms f t , with f 0 = Id, and has the form, where ∆ t and ϕ j (t) are the pullbacks of ∆ t and ϕ j (t) to Ω 0 , respectively. Taking the variation gives˙ ∆ϕ j (0) + (∆ 0 + λ 2 j )˙ ϕ j (0) = 0. Take the inner product with ϕ k (0) in the same eigenspace. Integration by parts in the second term kills the second term. Thus we get ∆ ϕ j (0), ϕ k (0) = 0.
The variation∆ can be calculated (see for example [HeZe12]) to obtain: ∂Ω 0ρ (∂ ν ϕ j )(∂ ν ϕ k )ds = 0, for all ϕ j , ϕ k in the λ j -eigenspace of the Dirichlet problem. A similar proof works for the relevant quadratic form for the Neumann problem.

Proof of Theorem 4
Before we prove our main theorem, we need to study the limits of the equations (41) along sequences of eigenvalues introduced in Theorem 1.
We recall that η is the symplectic dual of the arclength variable s. From the equation ηds = p ϑ dϑ, we find that in the (ϑ, p ϑ ) coordinates, η is given by η = p ϑ c 2 (cosh 2 ρ max − cos 2 ϑ) .

SPECTRAL RIGIDITY OF THE ELLIPSE 19
The corollary follows in the Neumann case by taking out the constant α − cosh 2 ρ max from the integral (42).
Proof using invariant curves inside the separatrix. We change variables to u = cos ϑ and also set x = √ α. Then the integral (43) becomes Writing f (u) = P (u 2 ) √ 1−u 2 , this becomes

The transform
Af (x) = x 0 f (u) √ x 2 − u 2 du is closely related to the Abel transform. We claim that the left inverse Abel transform is given by, The key point is the integral identity, It follows that if Bg(u) is the integral in the purported inversion formula, u 0 x 0 Since A is left invertible, it follows that ker A = {0}. Since f (u) = P (u 2 ) √ 1−u 2 lies in its kernel, we have P = 0 and henceρ = 0.
Proof using invariant curves outside the separatrix. The proof of the second assertion of Proposition 19 is similar to the final steps in the proofs of spectral rigidity results of [GuMe79], [HeZe12], and [Vi20], for the ellipse in various settings. We need to show that (44) implies P = 0. We change variables by u = cos 2 ϑ and this time we set f (u) = P (u) √ u(1−u) . Then Since the left hand side as a function of α is smooth at cosh 2 ρ max , all its Taylor coefficients at this point must vanish. Thus By the Stone-Weierstrass theorem, f = 0, hence P = 0.
6.1. Infinitesimal rigidity and flatness. In Section 3.2 of our earlier paper [HeZe12], we proved that infinitesimal rigidity implies flatness, which completes the proof of Corollary 5: