{'fraction_non_alphanumeric': 0.05066442859141139, 'fraction_numerical': 0.013372499650300742, 'mean_word_length': 3.982992960200739, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spectral rigidity for spherically symmetric manifolds with boundary

We prove a trace formula for three-dimensional spherically symmetric Riemannian manifolds with boundary which satisfy the Herglotz condition: The wave trace is singular precisely at the length spectrum of periodic broken rays. In particular, the Neumann spectrum of the Laplace--Beltrami operator uniquely determines the length spectrum. The trace formula also applies for the toroidal modes of the free oscillations in the earth. We then prove that the length spectrum is rigid: Deformations preserving the length spectrum and spherical symmetry are necessarily trivial in any dimension, provided the Herglotz condition and a generic geometrical condition are satisfied. Combining the two results shows that the Neumann spectrum of the Laplace--Beltrami operator is rigid in this class of manifolds with boundary.


Introduction
We establish spectral rigidity for spherically symmetric manifolds with boundary. We study the recovery of a (radially symmetric Riemannian) metric or wave speed rather than an obstacle. To our knowledge it is the first such result pertaining to a manifold with boundary. We require the so-called Herglotz condition while allowing an unsigned curvature; that is, curvature can be everywhere positive or it can change sign, and we allow for conjugate points. Spherically symmetric manifolds with boundary are models for planets, the preliminary reference Earth model (PREM) being the prime example. Specifically, restricting to toroidal modes, our spectral rigidity result determines the shear wave speed of Earth's mantle in the rigidity sense.
The method of proof relies on a trace formula, relating the spectrum of the manifold with boundary to its length spectrum, and the injectivity of the periodic broken ray transform. Specifically, our manifold is the Euclidean annulus M =B(0, 1) \B(0, R) ⊂ R n , R > 0 and n ≥ 2, with the metric g(x) = c −2 (|x|)e(x), where e is the standard Euclidean metric and c : (R, 1] → (0, ∞) is a function satisfying suitable conditions. (In appendix C we clarify that any spherically symmetric manifold is in fact of the form we consider -radially conformally Euclidean.) Our assumption is that the Herglotz condition d dr (r/c(r)) > 0 is satisfied everywhere. This condition was first discovered by Herglotz [11] and used by Wiechert and Zoeppritz [18].
By a maximal geodesic we mean a unit speed geodesic on the Riemannian manifold (M, g) with endpoints at the outer boundary ∂M := ∂B(0, 1). A broken ray or a billiard trajectory is a concatenation of maximal geodesics satisfying the reflection condition of geometrical optics at both inner and outer boundaries of M . If the initial and final points of a broken ray coincide at the boundary, we call it a periodic broken ray -in general, we would have to require the reflection condition at the endpoints as well, but in the assumed spherical symmetry it is automatic. We will describe later (definition 4.1) what will be called the countable conjugacy condition which ensures that up to rotation only countably many maximal geodesics have conjugate endpoints. This is a "generic" condition for reasonable wave speeds.
The length spectrum of a manifold M with boundary is the set of lengths of all periodic broken rays on M . If M is a spherically symmetric manifold as described above, we may choose whether or not we include the rays that reflect on the inner boundary r = R. If the radial sound speed is c, we denote the length spectrum without these interior reflections by lsp(c) and the one with these reflections by lsp (c). See figure 1 below for an illustration of these two types of geodesics. If the inner radius is zero (R = 0), the manifold is essentially a ball and the two kinds of length spectra coincide. We note that every broken ray is contained in a unique two-dimensional plane in R n due to symmetry considerations. Therefore, it will suffice to consider the case n = 2; the results regarding geodesics and the length spectrum carry over to higher dimensions.
We denote the Neumann spectrum of the Laplace-Beltrami operator in three dimensions, ∆ c = c 3 ∇ · c −1 ∇, on M by spec(c), where we impose Neumann-type boundary conditions on both the inner and outer boundary. The spectrum spec(c) includes multiplicity, not just the set spectrum.
The result holds true also for the length spectrum lsp (c) including reflections at the inner boundary, provided that R > 0. Theorem 1.2. Let B =B(0, 1) \B(0, R) ⊂ R 3 , R ≥ 0, be an annulus (or a ball if R = 0). Fix ε > 0 and let c τ , τ ∈ (−ε, ε), be a C ∞ function [R, 1] → (0, ∞) satisfying the Herglotz condition and the countable conjugacy condition and depending C 1 -smoothly on the parameter τ . Assume also that the length spectrum is non-degenerate in the sense that if two periodic broken rays have the same primitive period, then they differ by a rotation. If R = 0, we assume that all odd order derivatives of c τ vanish at 0. If spec(c τ ) = spec(c 0 ) for all τ ∈ (−ε, ε), then c τ = c 0 for all τ ∈ (−ε, ε). Theorem 1.1 is restated as theorems 4.6 (for lsp(c)) and 4.7 (for lsp (c)). We also prove an analogous theorem for the union of various spectra of this kind (theorem 4.8). We note that the dimension is irrelevant for the length spectral rigidity results; if the sound speed is fixed, the length spectrum is independent of dimension.
Using proposition C.1, we find the following corollaries: Corollary 1.3. Let B =B(0, 1) \B(0, R) ⊂ R n , n ≥ 2 and R ≥ 0, be an annulus (or a ball if R = 0). Fix ε > 0 and let g τ , τ ∈ (−ε, ε), be a C 1,1 -regular non-trapping SO(n)-invariant Riemannian metric making the boundary strictly convex and satisfying the Figure 1. The first figure shows a closed orbit whose period is in lsp(c) and lsp (c). The second figure shows a closed orbit whose period is in lsp (c), but not in lsp(c).
In dimension n ≥ 3 all rotation equivariant diffeomorphisms (diffeomorphisms commuting with the action of SO(n)) are radial. In dimension n = 2 the diffeomorphisms are also radial if the metrics are assumed to be O(2)-invariant. Corollary 1.4. Let B =B(0, 1) \B(0, R) ⊂ R 3 , and R ≥ 0, be an annulus (or a ball if R = 0). Fix ε > 0 and let g τ , τ ∈ (−ε, ε), be a C ∞ -regular non-trapping rotation invariant Riemannian metric making the boundary strictly convex and satisfying the countable conjugacy condition and depending C 1 -smoothly on the parameter τ . Assume also that the length spectrum is non-degenerate up to rotations. If spec(g τ ) = spec(g 0 ) for all τ ∈ (−ε, ε), then there are radial diffeomorphisms ψ τ : B → B so that ψ * τ g τ = g 0 for all τ ∈ (−ε, ε). Remark 1.5. While our main results are stated for the Laplace-Beltrami operator, they are equally valid for the spectra associated to the toroidal modes (see [6, chapter 8.6] for a precise definition) of an elastic operator. In section 2, we point out that the length spectrum may be recovered from either the spectrum of the Laplace-Beltrami operator or the toroidal mode eigenfrequencies so the results above hold when one has the spectrum associated to toroidal modes. In fact, our trace formula (proposition 2.2) to recover the length spectrum from the spectrum is done for the more general toroidal modes and frequencies.
Nevertheless, we show how our proof and result holds for the Laplace-Beltrami operator as well. We note that the spectrum of the Laplace-Beltrami operator depends on dimension.
The proofs of these theorems require three key ingredients which we elaborate in the next subsection: • A trace formula relating the Neumann spectrum of the Laplacian to the length spectrum. Thus, we will only need to prove length spectral rigidity since spectral rigidity (theorem 1.2) follows from the trace formula. • Sufficiently many periodic broken rays are stable under geometry-preserving perturbations of the metric and the derivative of the length of such broken rays is the periodic broken ray transform of the variation of the metric. • The periodic broken ray transform uniquely determines a radial function.
Outline of the proof. The breakup of the paper is as follows. In section 2, we discuss the relevant partial differential operators and their eigenfunctions. We also discuss geodesics in spherical symmetry and state the trace formula that we will prove (see proposition 2.2). Section 3 will be devoted to a proof of the trace formula. It will be necessary to use several transforms and the Debye expansion (see Appendix A) to convert the Green's function for the wave propagator written in terms of eigenfunctions to the dynamical Fourier integral operator (FIO) representation analogous to the one in [8] and [10]. Along the way, the connection from eigenfunction to geodesics becomes rather explicit. To reinforce this point, in section 3.3 we revisit the wave propagator constructed in [10] and show how all our explicit calculations relate to the abstract, geometric construction in that paper. Finally, a standard application of the method of steepest descent and stationary phase provides the leading order asymptotics for the trace.
After proving the trace formula, we prove the rigidity of the length spectrum in section 4. Together with the trace formula this implies the rigidity of the Neumann spectrum of the Laplace-Beltrami operator. We have a family of radially symmetric wave speeds c τ parameterized by τ ∈ (−ε, ε). For any periodic broken ray that is locally stable under the family of deformations, the derivative of its length is the integral of the metric variation over the periodic broken ray. In the case of closed manifolds and periodic geodesics this is well known, and in the case of non-periodic broken rays this was observed in [12]. Since the length spectrum is independent of the parameter τ , these derivatives vanish. The countable conjugacy condition (definition 4.1) guarantees that sufficiently many periodic broken rays are stable, so that we may conclude that the periodic broken ray transform of the function d dτ c −2 τ τ =0 vanishes. It then follows from recent results of periodic broken ray tomography on spherically symmetric manifolds [7] that the function in question has to vanish. Since the function is radial, this can be seen as an injectivity result for an Abel-type integral transform. Consequently c τ is independent of τ and the rigidity of the length spectrum and thus the spectrum follows.
Acknowledgements. M.V.d.H. gratefully acknowledges support from the Simons Foundation under the MATH + X program and the National Science Foundation under grant DMS-1559587. J.I. was supported by the Academy of Finland (decision 295853), and he is grateful for hospitality and support offered by Rice University during visits. V.K. thanks the Simons Foundation for financial support. We would also like to thank Gunther Uhlmann for helpful discussions and providing us useful references.

Geodesics and Eigenfunctions
In this section, we describe the relevant partial differential operators, associated eigenfunctions, and the connection between toroidal modes and the spectrum of the Laplace-Beltrami operator. We state our trace formula (proposition 2.2) along with an important remark related to the Laplace-Beltrami operator described in the introduction. First, let us provide a preliminary discussion of geodesics in spherically symmetric manifolds.
2.1. Geodesics in a spherically symmetric model. For the moment, we suppose that n = 2 and equip the annulus M =B(0, 1) \B(0, R) with spherical coordinates. For a maximal geodesic we define its radius as its Euclidean distance to the origin. We let γ(r) be the maximal geodesic of radius r which has its tip (the closest point to the origin) at the angular position θ = 0.
For r 0 ∈ (R, 1), the geodesic γ(r 0 ) can be parametrized as so that r(0) = r 0 and θ(0) = 0. Here L(r 0 ) > 0 is the half length of the geodesic. Using the conserved quantities c(r(t)) −2 [r (t) 2 + r(t) 2 θ (t) 2 ] = 1 (squared speed) and p = p γ := c(r(t)) −2 r(t) 2 θ (t) = r 0 /c(r 0 ) (angular momentum) one can find the functions r(t) and θ(t) explicitly. Using these conserved quantities it is straightforward to show that where β(r; p) := c −2 (r) − r −2 p 2 . We introduce the quantity β because it will appear naturally in the asymptotic approximations of eigenfunctions in section 3 and its relation to geodesics will now be clear. We denote α(r) := θ(L(r)), where θ is the angular coordinate (taking values in R) of the geodesic γ(r). That is, 2α(r) is the angular distance of the endpoints of γ(r). It may happen that α(r) > π if the geodesic winds around the origin several times. Using the invariants given above, one can also find an explicit formula for α(r): dr .
We will use the following lemma without mention whenever we need regularity of these functions.
Lemma 2.1. When c is C 1,1 and satisfies the Herglotz condition, then the functions α and L are C 1 on (R, 1).
2.2. Toroidal modes, eigenfrequencies, and traces. We now use spherical coordinates (r, θ, ψ). Toroidal modes are precisely the eigenfunctions of the isotropic elastic operator that are sensitive to only the shear wave speed. We forgo writing down the full elastic equation, and merely write down these special eigenfunctions connected to the shear wave speed (full details with the elastic operator may be found in [6, chapter 8.6]). Analytically, these eigenfunctions admit a separation in radial functions and real-valued spherical harmonics, that is, in which k = l(l + 1) (instead of the asymptotic Jeans relation, k = l + 1 2 ) and U represents a radial function ( n U l ). In the further analysis, we ignore the curl (which signifies a polarization); that is, we think of n D l as the multiplication with n U l (−k −1 ). In the above, Y m l are spherical harmonics, defined by The function, U (a component of displacement), satisfies the equation where µ = µ(r) is a Lamé parameter and ρ = ρ(r) is the density, both of which are smooth. Also, ω = n ω l denotes the associated eigenvalue. Here, l is referred to as the angular order and m as the azimuthal order. The traction is given by which vanishes at the boundaries (Neumann condition). The radial equations do not depend on m and, hence, every eigenfrequency is degenerate with an associated (2l+1)-dimensional eigenspace spanned by {Y −l l , . . . , Y l l }. We use spherical coordinates (r 0 , θ 0 , ψ 0 ) for the location, x 0 , of a source, and introduce the shorthand notation ( n D l ) 0 for the operator expressed in coordinates (r 0 , θ 0 , ψ 0 ). We now write the (toroidal contributions to the) fundamental solution as a normal mode summation On the diagonal, (r, θ, ψ) = (r 0 , θ 0 , ψ 0 ) and, hence, Θ = 0. Here Θ is the angular epicentral distance, cf. (3.2). We observe the following reductions in the evaluation of the trace of (2.5): • The functions U (r) are normalized, so that Meanwhile, the spherical harmonic terms satisfy l m=−l Y m l (θ, ψ) 2 sin θ dθ dψ = 2l + 1 (counting the degeneracies of eigenfrequencies). • If we were to include the curl in our analyis (generating vector spherical harmonics), taking the trace of the matrix on the diagonal yields sin θ dθ dψ = 2l + 1.
From the reductions above, we obtain We now write n f l (t) = Re e i nωlt i n ω l which is the inverse Fourier transform of Moreover, taking the Laplace-Fourier transform yields This confirms that the trace is equal to the inverse Fourier transform of ∞ l=0 (2l + 1) ∞ n=0 1 2 i n ω l πδ(ω − n ω l ) − πδ(ω + n ω l ) .
In the next subsection, we explain how the toroidal eigenfrequencies { n ω l } n,l relate to the Neumann spectrum of the Laplace-Beltrami operator described in the introduction. We also show why all our results and proofs in connection to the trace formula (proposition 2.2) hold for this spectrum as well.
2.3. Connection between toroidal eigenfrequencies, spectrum of the Laplace-Beltrami operator, and the Schrödinger equation. We relate the spectrum of a scalar Laplacian, the eigenvalues associated to the vector valued toroidal modes, and the trace distribution ∞ l=0 ∞ n=0 (2l + 1) cos(t n ω l ). We note that (2.3) and (2.4) is a 0th order operator and ω 2 is a particular eigenvalue. Hence U Y m l are scalar eigenfunctions for the self-adjoint (with respect to the measure ρ dx) scalar operator P with Neumann boundary conditions (on both boundaries) expressed in terms of N .
The above argument shows that we may view the toroidal spectrum { n ω 2 l } n,l as also the collection of eigenvalues λ for the boundary problem on scalar functions (2.8). Thus (2.2) can be written in the form where the last sum is taken with multiplicities for the eigenvalues. (While G is a vector valued distribution, the asymptotic trace formula we obtain is for Tr(∂ t G), which is equal to λ∈spec(P ) cos(t √ λ) by the normalizations we have chosen.) Up to principal symbols, P coincides with the ∆ c = c 3 ∇ · c −1 ∇ upon identifying c 2 with ρ −1 µ. This means that the length spectra of P and ∆ c will be the same even though they have differing subprincipal symbols and spectra. Thus, the trace formula which will appear to have a unified form, connects two different spectra to a common length spectrum and the proof is identical for both.
For concreteness, we recall [10, theorem 1]. Suppose that λ 1 ≤ λ 2 ≤ · · · → ∞ denote the Neumann spectrum of the Laplace-Beltrami operator ∆ g . We form the distribution (2.9) λ∈spec(∆g) Under certain geometric conditions (where there is no symmetry) for a simpler manifold described there, the authors prove the following: Let T be the singular support of (2.9) and the only closed geodesics, γ, of period T satisfy certain geometric conditions and have Maslov indices σ γ . Then for t near T , (2.9) is equal to the real part of (2.10) where T is the primitive period and P γ is a certain Poincaré map described there. Our results will settle the question of whether there is a formula analogous to (2.10) (same as [10, (1.3)]) for the distributions λ∈spec(P ) cos(t √ λ) and λ∈spec(c) cos(t √ λ) in our spherically symmetric manifold with boundary, encompassing a ball and an annulus 1 .
We will prove a trace formula using a WKB expansion of eigenfunctions. To this end, it is convenient to establish a connection with the Schrödinger equation. Indeed, we present an asymptotic transformation finding this connection. In boundary normal coordinates (r, θ) (which are spherical coordinates in dimension three by treating θ as coordinates on the 2-sphere), where ∆ θ is the Laplacian on the 2-sphere. We write λ = ω 2 as before and rewrite the second-order equation (2.8) The Neumann condition is applied at r = 1 and at r = R. In preparation of the asymptotic analysis we will instead invoke the transformation using the same notation (by abuse of notation, cf. (2.11)) With this definition, the matrix A in (2.12) admits the expansion viewing ω as a large parameter. A similarity transform gives We then seek an asymptotic transformation L and W so that the original matrix system implies ∂ r W = ωΛ(r)W.
Substitution and equating terms with equal powers of ω −1 gives and so on. We find a simple solution Then Here, w satisfies the equation If y(θ) is an eigenfunction of ∆ θ with eigenvalue −k 2 and V = V (r) is radial function, we choose w = w k = V (r)y(θ) so that V (r) must now satisfy generating two linearly independent solutions. The WKB asymptotic solution to this PDE with Neumann boundary conditions will precisely give us the leading order asymptotics for the trace formula, and is all that is needed. We note that rU = µ −1/2 V and rT (r) = µ 1/2 ∂ r V.
For the boundary condition, we note that we would end up with the same partial differential equation with different boundary conditions for V in the previous section if we had used the boundary condition ∂ r u = 0 on ∂M . Indeed, one would merely choose N u = µ∂ r u instead without the 0'th order term. However, the boundary condition for V would be of the form with K signifying a smooth radial function. Nevertheless, the leading order (in ω) asymptotic behavior for V stays the same despite the K term as clearly seen in the calculation of section 3.1.1. Thus, our analysis applies with no change using the standard Neumann boundary conditions. This should come as no surprise since in [10], the 0'th order term in the Neumann condition played no role in the leading asymptotic analysis of their trace formula. Only if one desires the lower-order terms in the trace formula would it play a role.

2.4.
Poincaré maps and the trace formula. Here, we describe the relevant Poincaré map that will appear in the trace formula and state the trace formula that we will prove. Let R t → γ(t) be a periodic broken bicharacteristic in S * M of period T > 0 (see [10] for the relevant definitions). It is associated to the metric c −2 (|x|)e, where c is a smooth radial function, e is the Euclidean metric, and γ undergoes reflections in ∂S * M according to Snell's law. We also denote by Φ T : S * M → S * M the broken bicharacteristic flow of T units of time as described in [10]. Its fixed point set is given by and without loss of generality, we assume that C T is connected, for otherwise, we would look at a connected component instead. We impose the clean intersection hypothesis appearing in [8,10] so that C T is a submanifold for any T and T m C T = ker(Id −dΦ T (m)) at each point m ∈ C T . This holds for all periodic orbits if and only if c satisfies the periodic conjugacy condition (definition 4.2) which requires that the endpoints of a single maximal geodesic segment of a periodic broken ray are never conjugate; see remark 4.3. By construction, the image of γ belongs to C T . There is an obvious symplectic group action of SO(3) on T * M under which the Hamiltonian 1 2 c 2 |ξ| 2 is invariant; here, ξ denotes the dual variable to x. Thus, for each g ∈ SO(3), the set g · Im(γ) given by the group action also belongs to C T . Assuming that c has no other symmetries (follows from the Herglotz condition) and the periodic conjugacy condition, all elements of C T are obtained this way. This is because a periodic orbit will either fail to be periodic or not have period T after a small perturbation of the angular momentum p; hence, p remains constant on C T . Thus, C T may be parameterized by R t × SO (3), revealing that under the Herglotz and periodic conjugacy conditions The Herglotz condition ensures that the group action never coincides with the geodesic flow; without it, the dimension could quite possibly be smaller. Hence, T m S * M/T m C T is only one-dimensional and we obtain an induced map on the quotient space, We denote the equivalence class of all closed orbits of period T related to γ by an element of SO(3) or by a time reversal of γ by [γ]. We write the above map as I − P [γ] and refer to P [γ] as the Poincaré map associated to the equivalence class of γ. Now, I − P [γ] will end up being an isomorphism, whose determinant at each point m ∈ C T will stay invariant. Hence, the quantity I − P In the above, γ may be multiple revolutions of another closed orbit of minimal period called the primitive orbit associated to γ, which has a primitive period denoted T . Note that T will merely be a positive integer multiple of T . In spherical symmetry, γ is confined to a disk and it must be a concatenation of geodesic segments that travel from the outer boundary r = 1 to either the reflection point r = R or the turning point r = R (see section 4 for details). We let N γ denote the number of these segments comprising the primitive orbit associate to γ.
We now state our proposition pertaining to the trace formula.
Proposition 2.2. Suppose the radial wave speed c satisfies the Herglotz condition and the periodic conjugacy condition (definition 4.2). Assume additionally that the length spectrum is non-degenerate in the sense that any two periodic broken rays of the same length are rotations of each other. The distribution (Tr (∂ t G))(t) = n,l (2l + 1) cos(t n ω l ) is singular precisely at the length spectrum. Suppose T ∈ singsupp(Tr(∂ t G)) and let d be the dimension of the fixed point set for Φ T . Suppose that γ is one of the broken periodic orbits of period T . Then d = 4 and for t near T , the contribution of [γ] to the leading singularity of (Tr(∂ t G))(t) is the real part of In appendix B we identify this trace formula in the framework of manifolds with symmetries given by a compact Lie group. In appendix D we discuss some edge cases that justify our geometric assumptions. Remark 2.3. This trace formula is in fact a more general statement than that for the Neumann Laplace-Beltrami operator, which is just a special case of the above proposition.
Remark 2.4. The periodic conjugacy condition and non-degeneracy of the length spectrum are needed to prove that the singular support of the trace is precisely the length spectrum. However, they are not necessary for spectral rigidity or proving that the spectrum determines the length spectrum.
For unique determination of the length spectrum, it is enough that the primitive length spectrum (excluding all but primitive orbits) is non-degenerate. Given the singularity at T , we know what the singularities at 2T , 3T , . . . will be. If they are not as expected, then another broken ray must contribute a singularity at the same place, and we have found the next primitive length. This allows to recover the primitive length spectrum and therefore the whole length spectrum from the (shapes and locations of) singularities in the trace. But if two primitive lengths coincide, there is no way to distinguish the corresponding singularities.
If we drop the periodic conjugacy condition, then some periodic orbits may fail the clean intersection hypothesis. This allows us to recover only a part of the length spectrum from the singularities, but this part is enough. Such problematic periodic broken rays are ignored anyway in the proof of length spectral rigidity, since they might not be stable under deformations.

Proof of the Trace Formula (proposition 2.2)
In this section, we prove the trace formula in the form of proposition 2.2 for the annulus. The idea behind the proof is to construct rather explicitly a fundamental solution. First, we do some preliminary analysis to manipulate G into the right form before taking its trace. Concretely, in subsection 3.1, we construct WKB eigenfunction solutions to get explicit formulas for the leading order asymptotics of the eigenfunctions. Afterwards, in subsection 3.2 we use the classical Poisson summation formula and the Debye expansion to write the leading order asymptotics for G as a certain propagator, which relates the eigenfunctions to geodesics in the annulus. At that point, we use section 3.3 to show how all our constructions are quite natural and directly relate to the wave propagator appearing in [10]. Finally, we complete the proof in section 3.4 by taking traces and carrying out the method of steepest descent and stationary phase to obtain the desired asymptotic formula appearing in proposition 2.2.
In the further analysis, we employ the summation formula, where the P l are the Legendre polynomials, with P l (1) = 1, and Θ signifies the angular epicentral distance, Remark 3.1. We note thatĜ is the kernel of the resolvent in the time-harmonic formulation. The normal mode summation becomes explicitly showing the eigenfrequencies as simple poles (cf. (2.7)).
We abuse notation and denote in the formula for G to not treat the curl operations at first. This will not cause a risk of confusion since we will specify the exact moment that we apply the curl operators, which will be just before taking the trace in subsection 3.4.
3.1. Asymptotic analysis of eigenfunctions. We describe the radial eigenfunctions and their asymptotic expansions for general k ∈ R. We then introduce the dispersion relations, ω n (k).
3.1.1. WKB eigenfunctions. Here, we consider asymptotic solutions, V , to (2.15). Depending on p, we distinguish the following regimes: Here, β 2 (r) < 0, and the solution is always non-oscillatory, that is, evanescent. We do not obtain eigenfunctions. • Diving (R/c(R) < p < 1/c(1)): We summarize the WKB solution of (2.15) in the vicinity of a general turning point. A turning point, r = R , is determined by Near a turning point, r ≈ R , and Away from a turning point, Matching asymptotic solutions yields From these one can obtain a uniform expansion, that is, the Langer approximation V (r, ω; p) = 2π 1/2 χ 1/6 (−β 2 ) −1/4 Ai(χ 2/3 (r)), valid for r ∈ [R, 1]. One obtains eigenfunctions corresponding with turning rays. Up to leading order, where r R , Here, B is obtained from the normalization (2.6), which requires the uniformly asymptotic solution over the entire interval [R, 1]. Applying the Riemann-Lebesgue lemma, one obtains Here, The solutions are oscillatory in the entire interval [R, 1] (β 2 (r; p) > 0), correspond with reflecting rays, and are of the form to leading order. Imposing the Neumann condition, T (R) = 0, implies that D = C. The constant C is obtained from the normalization (2.6) using the oscillatory solution over the entire interval [R, 1]. Applying the Riemann-Lebesgue lemma yields Here, 1 R β −1 c −2 dr can be identified with the half one-return reflection travel time, 3.1.2. Boundary condition and dispersion relations. We backsubstitute p = ω −1 k in β. Imposing the Neumann boundary condition yields: • Diving (R/c(R) < p < 1/c(1)): π (n is the overtone index).
All these (radial) quantization-type conditions yield solutions n ω k =: ω n (k). Using the implicit function theorem, we introduce k n = k n (ω) as the solution of We revisit the general relation between phase and group velocities. We have The corresponding ray parameter is given by
3.2.1. Analytic continuation. The Legendre equation is of second order and hence admits two linearly independent solutions: Legendre functions of the first and second kind, P λ , Q λ , with integral representations 1]. These integral representations are valid for λ > −1. Analytic continuation of the integral representation for P λ , Q λ into the region Re{λ} ≤ −1 follows the relations The bottom equation implies that Q λ has simple poles at the negative integers λ = −1, −2, . . .. The Legendre function of the first kind, P λ , coincides with the Legendre polynomials for λ = l = 0, 1, 2, . . . (cf. (3.1)).

Application of Poisson's formula.
Poisson's formula is given by Remark 3.2. Poisson's formula can be obtained from the Watson transformation: If f is a function analytic in the vicinity of the real axis, and C = C − ∪ C + is a contour around the positive real axis, then The integrand in the right-hand side has simple poles at (2n + 1)/2, n ∈ N 0 -where cos(kπ) = 0. It follows from the residual theorem. Poisson's formula is obtained as follows.
In the limit s → ∞ the path of integration is C − , while in the limit s → −∞ the path of integration is C + . One then expands [cos(kπ)] −1 = 2e − i kπ 1 + e −2 i kπ in a series separately for Im{k} > 0 and < 0.
We apply Poisson's formula to the summation in l in (3.3) while keeping the summation in n intact. We obtain This expression can be rewritten aŝ Traveling-wave Legendre functions. Traveling-wave Legendre functions are given by Using (3.9)-(3.10) yields the continuation relations Both have simple poles at the negative integers. Remark 3.3. Their asymptotic behaviors, as |λ| 1, are (assuming that Θ is sufficiently far away from the endpoints of [0, π]) upon substituting z = cos Θ. Taking into consideration the time-harmonic factor e i ωt , it follows that Q (1) represents waves traveling in the direction of increasing Θ, while Q (2) represents waves traveling in the direction of decreasing Θ.
Substituting (3.17) into (3.13), we get In the application of (3.11)-(3.12), the contour C − is followed for the series containing e − i(s−1)kπ and e − i skπ , while the contour C + is followed for the series containing e i(s+1)kπ and e i(s−2)kπ .

3.2.3.
Preparation for the method of steepest descent. In preparation of the application of the method of steepest descent, we rewrite the k-integral from R ≥0 to R. We make use of the analytic extension of Q k−1/2 (cos Θ) + e ± i kπ tan(kπ)P k−1/2 (− cos Θ).
Using the symmetry off n (.; ω), H n (.), it follows that the integrals over P k−1/2 (cos Θ) cancel. Then The integrands in the terms of these series can be identified as wave constituents travelling along the surface or boundary, the representations of which can be obtained by techniques from semi-classical analysis. Indeed, s is referred to as the (multi-orbit) arrival number while we distinguish the orientation of propagation in the two series. The term s = 1 corresponds with waves that propagate from source to receiver along the minor arc; the term s = 2 corresponds with waves that propagate from source to receiver along the major arc. At Θ = 0 and Θ = π the traveling wave Legendre functions have logarithmic singularities, namely log Θ at Θ = 0 and log(π − Θ) at Θ = π; the singularities cancel when taking the sums together.

3.2.4.
Traveling wave expansion. Here, we apply to (3.18) the Debye expansion described in appendix A to obtain a form more closely resembling a wave propagator: where the τ i are time-delay functions which we will relate to the travel time of a geodesic, and N i are integers contributing to the KMAH index as described in appendix A. Next, we change variables of integration from k to p. We encounter the Jacobian (cf. (3.8)) dp dk The path of integration is beneath the real axis, while taking ω > 0. After making the substitution into (3.18), we obtain ωp−1/2 (cos Θ)e − i ω(s−1)pπ p −1 dp ωp−1/2 (cos Θ)e − i ωspπ p −1 dp.
In the next section, we review the abstract argument appearing in [10] in order to justify why we consider the above expression a wave propagator. It will also motivate the method of steepest descent and stationary phase calculation that we perform in section 3.4.

3.3.
Melrose-Guillemin wave propagators. In [10], Guillemin and Melrose show how the Neumann half-wave propagator, e ± i t √ −∆g , may be written as a sum of operators denotedṼ i ± (t), such that for a fixed t, it is a canonical graph FIO whose canonical relation is a certain billiard map in phase space that we will briefly describe. In order to avoid corners, they embed M in a boundaryless manifoldM and thenṼ i + is an FIO onM × M . It maps a covector (y, η) ∈ S * M o to the "time t" endpoint of a broken bicharacteristic that undergoes i reflections at points x k (y, η) ∈ ∂M , k = 1, . . . , i. Precisely, let t i = t i (y, η) denote the travel time along the broken bicharacteristic to the last reflection point (x i , ξ i ) ∈ ∂ + S * M and let (x i , ξ r i ) ∈ ∂ − S * M be the reflected covector pointing "inside" M . The canonical graph maps (y, η) to Φ t−t i (x i , ξ r i ) ∈ S * M (here, Φ t denotes the bicharacteristic flow described in section 2). If fewer than i reflections take place by time t or it maps outside of M , thenṼ i + (t) is microlocally smoothing at such covectors. It will be convenient to denote a covector locally in polar coordinates as η = |η| g(y)η ∈ T * y M . Then the corresponding conic Lagrangian may by parameterized by a phase function of the form ϕ i (t, x, y, η) = |η| g (−t + S i (x, y,η)), where, essentially, S i gives the travel time between points x, y of the broken geodesic undergoing i reflections that starts at (y,η) (when one findsη ∈ S * y M that minimizes S i ). For the following calculations, all that is necessary is thatṼ i + (t) is a canonical graph FIO and that ϕ i is a phase function that locally parameterizes the conic Lagrangian associated to the canonical graph. Thus, the Schwartz kernel ofṼ i + (t) is indeed given locally (as described in [10]) by (2π) −n e i(−t|η| g +S i (x,y,η)) a(t, x, y, η) dη, where a is a classical symbol of order m = 0. If we introduced spherical coordinates in η with radial variable ω := |η| g and took the leading order (homogeneous of degree m = 0) part of the classical symbol a, then it becomes clear that (3.19) has the same form after a Fourier transform F ω→t with phase function −ω(τ i (r, r 0 ; p) + pΘ + (s − 1)pπ) corresponding to S i = ωS i (x, y,η) above. The order in ω does not match yet because we have not yet applied the curl operators nor ∂ t . The similarities become clearer as we proceed with the stationary phase calculations for both operators.
We change variables into polar coordinates by first defining a map f : Then f |dT * M | = ω n−1 dω |dS * M |. Without loss of generality, since we only consider the leading order terms, we may assume that a is an m-homogeneous (in η) symbol. After a change of variables into polar coordinates (keeping in mind that a is supported in a local coordinate chart where S * M is trivialized), the above integral is: x, x, p) dx dp.
Stationary phase with a degenerate phase function. To obtain the leading order asymptotics, we apply the method of stationary phase. We denote the critical set of S i (viewed as function on S * M ) locally as We may break up this set into a countable number of disjoint connected components given by the values of S i : By construction of the phase function, C T ik is precisely the fixed point set of the bicharacteristic flow Φ T ik . Let d ik be the codimension of C T ik . If C T ik intersects ∂M transversely, then we may find independent coordinates z = z 1 , . . . , z 2n−1 such that C T ik is given by the vanishing of z 1 , · · · , z d ik and ∂M is given by z 2n−1 = 0. For notation purposes, it will be convenient to denote z I = (z 1 , · · · , z d ik ), z II = (z d ik +1 , · · · , z 2n−2 ) and O a neighborhood where such coordinates are valid.
We assume that S i is nondegenerate in the directions normal to the critical manifold C T ik , so that after a change of coordinates still denoted by the same letters, S i is locally given by This follows precisely from the generalized Morse lemma. The Jacobian due to such a change of coordinates is identically 1 when restricted to C T ik and so we exclude it on our analysis to have simpler formulas.
In the above, we applied stationary phase in the variables z I corresponding to normal directions of C T ik . This requires the normal Hessian of S i which is non-degenerate: so that the leading amplitude after stationary phase becomes a (t, 0, z II , z 2n−1 ) = e i π 4 sgn(d 2 One may also show that J(0, z II , z 2n−1 ) dz II dz 2n−1 is the natural induced measure on C T ik obtained from dS * M in local coordinates. We label this measure as dµ ik . At this point, we may actually write Using integration by parts in ω, we may replace a(t, z) by a(T ik , z) since the difference will lead to a term that is an order lower in ω. Substituting this into the inner integral above, setting the leading order term in our trace formula becomes (with m = 0 since a is 0th order) In the case of isolated periodic geodesics, one has dim(C T ik ) = 1. The above calculation would then be the real part of For the spherically symmetric case, dim(C T ik ) = 4 which leads to a bigger singularity. In the next subsection, we apply analogous computations for taking the trace in our setting.

Proof of Proposition 2.2.
We are finally ready to apply the method of steepest descent and stationary phase to complete the proof of the trace formula.
Proof of proposition 2.2. We interchange the order of summation and integration, and invoke the method of steepest descent. We carry out the analysis for a single term, s = 1. For s = 2, 4, . . . we have to add spπ to τ i , and for s = 3, 5, . . . we have to add (s − 1)pπ to τ i , in the analysis below.
We find (one or more) saddle points for each i, where Later, we will consider the diagonal, setting r 0 = r and Θ = 0. We label the saddle points by k for each i (and s). We note that r, r 0 and Θ determine the possible values for p (given i) which corresponds with the number of rays connecting the receiver point with the source point (allowing conjugate points). For s = 1, the rays have not completed an orbit. With s = 3 we begin to include multiple orbits.
We carry out a contour deformation over the saddles and obtain The M ik contribute to the KMAH indices, while the T ik represent geodesic lengths or travel times. The orientation of the contour (after deformation) in the neighborhood of p k is determined by sgn ∂ 2 p τ i | p=p k . We note that • M ik = M s;ik (r, r 0 , Θ) for multi-orbit waves (s = 3, 4, . . .) includes polar phase shifts produced by any angular passages through Θ = 0 or Θ = π as well; • if r lies in a caustic, the asymptotic analysis needs to be adapted in the usual way.
Next, we apply both curl operations [−θ(sin θ) −1 ∂ ψ +ψ∂ θ ], [−θ 0 (sin θ 0 ) −1 ∂ ψ 0 +ψ 0 ∂ θ 0 ] to each term in the sum above and then tensor the vectors together in order to obtain a sum of 2-tensor fields. This will give us the actual normal mode summation of (3.3). Since we are interested in only the leading order asymptotics in ω, we need only consider these operations to the term exp[− i ωT ik ] which gives (ωp k ) 2 D(θ, ψ, θ 0 , ψ 0 ) where D is a 2-tensor field in the angular variables. We also apply an inverse Fourier transform followed by ∂ t to the formula above (since we are interested in the cosine propagator ∂ t G) to obtain to leading order It will now be convenient to denote the above quantity by F (1) if s is even and F (2) if s is odd. Thus, we have a sum of kernels of FIOs associated with the wave propagator. (Here, the summation over i, k signifies the summation over (broken) geodesics while the summation over s signifies the number of orbits, travelling clockwise or counterclockwise.) Since we will need to restrict ∂ t G to the diagonal (Θ = 0, r = r), we must be careful with the (sin Θ) −1/2 term. Fortunately, this term comes from the asymptotics of Q (1) ωp−1/2 and Q (2) ωp−1/2 which merely represent the two different directions of one particular geodesic. Of course, a periodic broken orbit has the same period no matter which direction one travels and in the trace formula, all orbits of a particular period are combined. Hence, we can combine Q (1) ωp−1/2 and Q (2) ωp−1/2 near Θ = 0, which is a logarithmic singularity that cancels when both of the functions are added together, and this will not affect the trace formula.
Precisely, after substituting the asymptotic formula for Q (1,2) ωp k −1/2 (cos Θ), we may write for s ≥ 3 odd Asymptotically, as Θ → 0 one has Q (1) ωp k −1/2 (cos(Θ)) 1 and so .)β(r 0 ; .)) −1/2 ∂ 2 p τ i (r, r 0 ; .) as Θ → 0. It will be convenient to denote A s;ik (r, r 0 , Θ) = (−) (s−1)/2 (rr 0 c(r)c(r 0 )) −1 (ρ(r)ρ(r 0 )) −1/2 · p(β(r; .)β(r 0 ; .)) −1/2 ∂ 2 p τ i (r, r 0 ; .) where now D is a scalar function, defined as the inner product of the two vector fields that make up D. One may check using l'Hospital's rule and equation (3.2) that D| θ=θ 0 ,ψ=ψ 0 = 2. Next, we take the trace of ∂ t G by restricting to (r = r 0 , Θ = 0) and integrating. The phase function on the diagonal is T ik = τ i (r, r, p k ) + π(s − 1)p k and we apply stationary phase in the variables r, θ, ψ with large parameter ω. Since one has ∂ p T i (r, r, p) = 0 at p = p k , the critical points occur precisely when After a quick calculation, the first condition forces T i to be independent of r. Also, we showed that for geodesics with turning points, U = O(ω −∞ ) when r < R . Finally, using the inverse Fourier transform, Setting R ik = {(r, θ, ψ); r ∈ [R, 1], ∂ r T i (r, r, p k ) + ∂ r 0 T i (r, r, p k ) = 0} whereR = R or R depending on p k , we have shown that in fact T ik remains constant over R ik so that only certain i are allowable. We find that modulo terms of lower order in ω, the trace microlocally corresponds with [13] Re is independent of r. We note that p k exists only for i, and s, sufficiently large, which reflects the geometrical quantization. From this expression, it is clear that the singular support of the trace consists of the travel times of periodic geodesics.

Remark 3.5.
It is now apparent how the above formula relates to the trace formula in [10]. A term above for a certain i, s, k index corresponds to the trace ofṼ ± i integrated over a critical manifold C T ik (which in our case is R ik ) as described in section 3.3. In both cases, the index i is used to keep track of the number of intersections of the broken ray with the boundary while the k index specifies a particular periodic ray and period.
The index s in T ik;s describing travel time has no analog for general manifolds and certainly does not appear in [8] and [10]. This is because the parameter s is merely a byproduct of using spherical coordinates to construct an FIO on the ball and using the particular phase function we have. It is used to keep track of the angular distance a particular geodesic has traversed in the disk, since, while the angular distance is greater than 2π when a particular geodesic traverses the full disk more than once, our angular coordinates only range from [0, 2π]. The factor (−) (s−1)/2 is part of the KMAH index corresponding to antipodal conjugate points in the sphere that the associated geodesic (when projected to r = 1) passes through. The other part of the KMAH index comes from M ik .
The Poincaré map appearing in [8] and [10] corresponds to the factor p −1 k ∂ 2 p τ i | p=p k −1/2 in our trace formula and the factor I − P [γ] −1/2 described in section 2. In [8], this factor generally varies over the entire critical manifold C T ik when its dimension is greater than 1, but will not in our case due to the symmetry. In fact, this is how we know that this factor corresponds to the Poincaré map: it must be a quantity that remains constant over the critical manifold R ik .
We further simplify the above formula, that is, the integral involving A s;ik . First, since T ik is independent of r, then so is τ i (r, r; p) = τ i (p). Thus, we may pull ∂ 2 p τ i out of the integral involving A s;ik precisely because we are integrating over a closed orbit: R ik A s;ik (r, r, 0)ρr 2 sin θ dr dθ dψ We recall that the travel time T for a piece of a geodesic from r 0 ro r is dr .
Hence, denoting T ik as the primitive period of the geodesic, we obtain where N ik is the number of geodesic segments fromR to 1 along the primitive orbit with length T ik , and |SO(3)| is the volume of SO(3) under a Haar measure [5]. Substituting these calculations, the leading order term in the trace formula is In the next section, we use the trace formula to prove our main spectral rigidity theorems stated in the introduction.

Proof of Spectral Rigidity
In this section we will prove that the length spectrum is rigid. By the trace formula this will imply that the spectrum of the Laplace-Beltrami operator is rigid. 4.1. Conjugacy conditions. The following condition will be convenient: Definition 4.1. We say that a C 1,1 sound speed c satisfies the countable conjugacy condition if there are only countably many radii r ∈ (R, 1) so that the endpoints of the corresponding maximal geodesic γ(r) are conjugate along that geodesic.
Assuming this condition will eventually imply that the length spectrum is countable. Throughout this paper "countable" includes empty and finite sets, but for the sake of brevity we shall not write "at most countable". The next condition is directly related to the clean intersection property discussed earlier. We will return to this condition in section 4.4.

Definition 4.2.
We say that the radial wave speed c satisfies the periodic conjugacy condition if α(r) ∈ πQ implies α (r) = 0. Restating geometrically, this means that if a broken ray is periodic, then the endpoints of a geodesic segment are not conjugate along the segment. If the second option is true, it will fail at each endpoint of the maximal interval on which α vanishes. Since lim r→1 α(r) = 0 when the boundary is strictly convex, such an endpoint exists. Therefore, assuming the Herglotz condition, the clean intersection property for all periodic broken rays is equivalent with the periodic conjugacy condition of definition 4.2. We point out that the function α is not well defined and the dimension of the fixed point set can be 3 if the Herglotz condition is violated.

Conditions for periodicity.
A geodesic can be extended into a broken ray. The geodesic segments of a broken ray are rotations of each other. It is easy to see that the broken ray corresponding to the geodesic γ(r) is periodic if and only if α(r) ∈ πQ. We want to understand the set of radii r ∈ (R, 1) for which this is the case.
Lemma 4.4. Let P ⊂ (R, 1) be the set of radii for which the corresponding broken ray is periodic. Let C ⊂ (R, 1) be the set of radii r for which the endpoints of the geodesic γ(r) are conjugate along γ(r).
• If C has empty interior, then P is dense.
Proof. The radius r ∈ (R, 1) parametrizes a family of geodesics. By rescaling the speed (from the originally assumed unit speed), we may assume that all geodesics are parametrized by [−1, 1]. Differentiating the geodesic γ(r) with respect to r gives a non-trivial Jacobi field, a variation of a geodesic γ(r). The value of the Jacobi field at the endpoints of the geodesic describes the movement of the endpoints of the geodesic under variation. On the other hand, α(r) gives the endpoint of the geodesic. Therefore the tangential component of the Jacobi field at the endpoint is α (r). It follows from the reparametrization that the component normal to the boundary vanishes. Thus if α (r) = 0, then the endpoints of γ(r) are conjugate along γ(r).
On the other hand, if the endpoints are conjugate, then there is a non-trivial Jacobi field vanishing at the endpoints. In dimension two there can only be a one-dimensional space of such Jacobi fields, so the Jacobi field in question must be symmetric in the time parameter of the geodesic. Then we can identify the Jacobi field as corresponding to a variation of the parameter r. Combining this with the previous observation shows that r ∈ C ⇐⇒ α (r) = 0.
This proves the first claim.
Let R < a < b < 1. To prove the second claim, it suffices to produce r ∈ (a, b) so that α(r) ∈ πQ; see the discussion right before the statement of this lemma. For a contradiction, assume that α(r) / ∈ πQ for all r ∈ (a, b). Since α is continuous, this implies that α is constant on (a, b). Thus α vanishes on (a, b), so C has interior -a contradiction.
Lemma 4.5. If the sound speed c ∈ C 1,1 satisfies the countable conjugacy condition (definition 4.1) and the Herglotz condition, then C and P are countable, C is closed and P is dense in (R, 1).
Proof. The countable conjugacy condition directly states that C is countable. Therefore C cannot have interior, and density of P follows from lemma 4.4. Since α is continuous, the preimage of zero under it, the set C, is closed.
It remains to show that P is countable. If it was uncountable, some level set of α would be uncountable. An uncountable set has uncountably many accumulation points and α vanishes at every accumulation point of a level set of α. This implies that C is uncountable, which is impossible.

4.3.
Length spectral rigidity. The length spectrum of a manifold M with boundary is the set of lengths of all periodic broken rays on M . If M is a spherically symmetric manifold as described above, we may choose whether or not we include the rays that reflect on the inner boundary r = R. If the radial sound speed is c, we denote the length spectrum without these interior reflections by lsp(c) and the one with these reflections by lsp (c). If the inner radius is zero (R = 0), the manifold is essentially a ball and the two kinds of length spectra coincide.
For clarity, we state the following three length spectral rigidity theorems separately.
Notice that dimension is irrelevant for the statements; if the sound speed is fixed, the length spectrum is independent of dimension.
The following theorem states that the same rigidity result is true for any finite disjoint union of manifolds of the types given in theorems 4.6 and 4.7.
Remark 4.9. It does not matter whether for every periodic broken ray only the primitive period is included in the length spectrum, or of its all integer multiples. The proofs of the three theorems above are the same in both cases.
Even more might be true, and we propose the following conjecture.
Conjecture 4.10. Under certain geometric hypotheses, a spherically symmetric manifold is uniquely determined by its length spectrum.
A verification of the conjecture would imply that such a manifold is uniquely determined by its spectrum under some geometric assumptions.
Remark 4.11. Theorem 4.8 may seem like an unnecessary generalization of the two preceding theorems, but it has geophysical significance. Consider a spherically symmetric model of the earth. It essentially consists of three different parts, an inner ball and two nested annuli. The full length spectrum of the Earth is the set of all periodic orbits for the different polarized waves. The statement in theorem 4.8 (with N = 2 and N = 1) is, however, incomplete in the sense that the coupling between different polarizations and transmission at boundaries are ignored. However, this is the best toy model for which rigidity is currently known.
Radial symmetry is an excellent approximation of the earth or planets in general. This symmetry is not exact, and unfortunately our method requires precise symmetry. Assuming radial symmetry is not merely a matter of technical convenience, but a truly necessary assumption. A key ingredient in the proof is that many periodic broken rays are stable under deformations of the metric. In spherical symmetry the broken rays are only stable under deformations that preserve the symmetry, otherwise they are typically unstable and the proof falls apart.
The countable conjugacy condition should be a generic property of sound speeds. It is a technical assumption we do not like to make, but it does not hinder the relevance for our planetary model. The Herglotz condition is crucial for the geometry of the problem. Without it the manifold would trap some geodesics inside and the geometry of broken rays would be very different. In the commonly used Preliminary Reference Earth Model (PREM) both pressure and shear wave speeds satisfy the Herglotz condition piecewise. Due to the layered structure of the Earth both wave speeds have jumps, whereas our result assumes C 1,1 regularity.

No inner reflections.
In this subsection we prove theorem 4.6. All the lemmas here are stated under the assumptions of the theorem. By decreasing ε slightly we can assume that c τ (1) is bounded away from zero and infinity without any loss of generality.
Let P τ ⊂ (R, 1) be the set of radii for which the corresponding broken ray -which is unique up to rotations -is periodic with respect to the sound speed c τ . A priori P τ depends on τ . If r ∈ P τ , the corresponding broken ray has n τ (r) reflections and winds around the origin m τ (r) times. We choose all periodic broken rays to have minimal period, so the natural numbers n τ (r) and m τ (r) are coprime. If α τ (r) is the angle we defined earlier (previously without the dependence on τ ), we have the identity πm τ (r) = n τ (r)α τ (r).
Proof. It follows from the Herglotz condition that r/c τ (r) ≤ 1/c τ (1), whence c τ (r) ≥ Rc τ (1). Since c τ (1) is uniformly bounded from below, the functions c τ are all uniformly bounded from below. We assumed the sound speeds to be uniformly bounded from above in the theorem, so for some constant C > 1 we have 1 C < c τ (r) < C for all τ and r. We write LHS RHS if there is a constant C independent of τ , r and s so that LHS < C · RHS. By LHS ≈ RHS we mean LHS RHS LHS.
It follows from the previous observation that Combining this with equations (2.2) and (2.1) gives the desired uniform estimates for α τ and L τ .
Lemma 4.13. If R = 0, then lim r→0 α τ (r) = π 2 and lim r→0 L τ (r) = 1 0 c τ (r) −1 dr. Proof. In the limit the corresponding maximal geodesic tends to a diameter of the ball. This limit is a radial geodesic and it is geometrically obvious that the angle and length corresponding to it are the limits stated above.
We denote by S τ the set of radii r ∈ (R, 1) for which α(r) ∈ πQ and α (r) = 0. Radii in this set will correspond to stable broken rays as we shall see in lemma 4.16 below.
Lemma 4.15. The set S τ is countable and dense in (R, 1).
Proof. Let C τ denote the set of zeros of α τ . (We used the same notation in lemma 4.4.) By lemma 4.5 we know that P τ is dense and countable and C τ is closed and countable. Therefore S τ = P τ \ C τ is dense and countable.
The next lemma shows that periodic broken rays are stable under variations of the radial sound speed. Periodic broken rays on a highly symmetric manifold are typically not stable under all variations of the metric, but we are only looking at variations that preserve the symmetry.
Lemma 4.17. Let ϕ : (−δ, δ) → (R, 1) for any δ > 0 be a C 1 function satisfying ϕ(τ ) ∈ P τ for all τ ∈ (−δ, δ). Let γ 0 : [0, T ] → M be a periodic broken ray with radius ϕ(0). Then Proof. A version of this lemma with non-periodic broken rays of finite length and arbitrary variations of the metric was given in [12, theorem 17]. (There is an error in the formula of the cited theorem; the factor 2 should be replaced with 1 2 .) Applying that result for conformal variations and broken rays that have the same starting and ending points gives the desired claim.
Lemma 4.18. If f : M → R is a continuous radially symmetric function (identified as a function f : (R, 1] → R) and r ∈ P 0 , then the integral of f over any periodic geodesic (with respect to sound speed c 0 ) of radius r is Proof. The proof is a simple calculation, similar to the one leading to equation (2.1).
Lemma 4.19. If Af (r) is the function of (4.1), then the map f → Af takes continuous functions to continuous functions and is injective on the space of continuous functions.
In the C ∞ setting relevant for the spectrum of the Laplace-Beltrami operator this was shown by Sharafutdinov [16]. With all these lemmas as ingredients, it is difficult not to prove theorem 4.6.
By lemma 4.17 this implies that the variation of the wave speed, f = d dτ c −2 τ τ =0 : M → R, integrates to zero over all periodic geodesics of radius r. This function f is radially symmetric, so we can think of it as a function (R, 1] → R. By lemma 4.18 we know that Equation (4.2) is true for a dense set of radii r ∈ (R, 1) by lemma 4.15, so it follows from lemma 4.19 that in fact f vanishes identically. We have found that d dτ c τ = 0 at τ = 0. The choice τ = 0 was in no way important to this argument, so in fact d dτ c τ = 0 for all τ ∈ (−ε, ε). This means that all sound speeds c τ indeed coincide.
A key step in the reasoning in the preceding proof can be stated as follows: A radially symmetric function is uniquely determined by its integrals over all periodic broken rays. There is even a reconstruction formula for this problem: [7,Remark 27]. Radial symmetry is important, as a general smooth function is not uniquely determined. Only the even part of the function can be recovered efficiently. This is in sharp contrast to the case of geodesic X-ray tomography on such manifolds where the ray transform has no kernel. See [7] for details.

4.3.2.
Inner reflections included. In this subsection we will prove theorems 4.7 and 4.8.
The only ingredient we need to add to the proof of theorem 4.6 is the following lemma. We will prove the lemma after showing how it completes the proofs of the theorems. Lemma 4.20. Let M =B(0, 1)\B(0, R) ⊂ R n , R > 1 and n ≥ 2, be an annulus. Equip M with a radially symmetric C 1,1 sound speed satisfying the Herglotz condition. Then the set of all lengths of periodic broken rays that reflect on the inner boundary r = R is countable.
We note that no assumption was made on conjugate points in the lemma.
Proof of theorem 4.7. We denote the set of all lengths of periodic broken rays with respect to sound speed c τ that reflect on the inner boundary by H τ . Then we have lsp (c τ ) = lsp(c τ ) ∪ H τ . As in the proof of theorem 4.6, we have a dense set of radii r ∈ (R, 1) for which there is a family of corresponding periodic broken rays varying continuously in τ ; this set was denoted by S 0 . Since lsp (c τ ) is independent of τ and each set H τ is countable, the lengths of the periodic broken rays in this family must in fact be independent of τ . The rest of the proof can be concluded as that of theorem 4.6. We only need the integrals of the variations of the sound speed over broken rays that do not hit the inner boundary.
Proof of theorem 4.8. We denote the individual Riemannian manifolds in question by M 1 , . . . , M N and M 1 , . . . , M N , equipped with their respective sound speeds. Each one of them has a countable length spectrum, so the length spectrum of the whole system is still countable. We can then use the argument presented in the proof of theorem 4.7 to conclude for each manifold separately that the sound speed has to be independent of the parameter τ .
Proof of lemma 4.20. Consider a geodesic joining the two boundaries; a periodic broken ray of the kind we need to study is a finite union of rotations and reflections of such a geodesic. We parametrize this geodesic with arc length starting from the inner boundary. Due to symmetry the geodesic is confined to a two-dimensional plane, and in this plane we may use polar coordinates. In these coordinates the geodesic is [0, T ] t → (r(t), θ(t)), with r(0) = R and r(T ) = 1. We denote ρ(r) = r/c(r).
Let the angle between the geodesic and the radially outward pointing normal vector be ω. (The metric is conformally Euclidean so we may use the Riemannian metric or the Euclidean metric to measure angles at a point.) We assume for the time being that ω ∈ (0, π 2 ). Since the geodesic has unit speed, we have cos ω = ρ(R)θ (0). With this information it is easy to see that the constant value of the angular momentum ρ(r(t)) 2 θ (t) is ρ(R) cos ω.
Using unit speed and the conservation of angular momentum one easily finds that the change in the angular coordinate over the geodesic is This angular difference depends on ω, but it is most convenient to think of it as a function of z := ρ(R) cos ω (the angular momentum). We denote this angular difference by β(z). Its geometrical meaning is the same as that of α(r) defined earlier, but we use a different letter to avoid confusion. Since an easy calculation gives This derivative is positive, so β is in fact a homeomorphism from [0, 1] to its image; the limits at 0 and 1 can be checked to be finite.
We may exclude radial geodesics. For non-radial geodesics the integral of (4.3) is nonsingular and differentiation under the integral sign is simple. For the similar result corresponding to the diving waves (the function α), the derivative is more complicated; cf. [7, proposition 15].
The broken ray corresponding to angular momentum z is periodic if and only if β(z) ∈ πQ. Since β is a homeomorphism between intervals, the set of angular momenta z corresponding to periodicity is countable. Therefore the set of corresponding lengths is also countable.

Spectral rigidity.
We are now ready to prove theorem 1.2. The proof is trivial at this point, but we record it explicitly for completeness.
Proof of theorem 1.2. Let us first consider the simplest case where c satisfies the periodic conjugacy condition.
The trace of the Green's function is determined by spec(c) through proposition 2.2. Since the trace as a function of t ∈ R is singular precisely at the length spectrum, spec(c) determines lsp (c). Then rigidity of the spectrum spec(c) follows from that of lsp (c); see theorem 4.7.
Let us then drop the periodic conjugacy condition. The Neumann spectrum spec(c) still determines the trace of the Green's function. As pointed out in remark 2.4, the singularities determine the part of the length spectrum corresponding to periodic broken rays satisfying the clean intersection property. Under the Herglotz condition and the countable conjugacy condition a periodic broken ray of radius r satisfies the clean intersection property if and only if α (r) = 0; see remark 4.3. In the notation of lemma 4.4 and denoting the length of a geodesic of radius r by (r), the problematic primitive lengths are precisely (r) for r ∈ C ∩ P . Since we assumed the length spectrum to be non-degenerate, none of the lengths (r) for r ∈ P \ C coincide with the problematic ones. Access to all radii r ∈ P \ C is sufficient for the proof of length spectral rigidity. Notice that C is exactly the same set that corresponds to possibly unstable periodic broken rays, and this data was ignored in the proof of length spectral rigidity anyway.

Appendix B. General framework for symmetries in a manifold
In this appendix, we consider more general situations for trace formulas when the manifold has symmetries given by a compact Lie group and we show how the quantities appearing in proposition 2.2 are special cases of a general framework. Our constructions here are inspired by the work in [4,5,2]. In particular, we provide general representation of the critical manifolds C T described in sections 2 and 3.3 under symmetry from a compact Lie group G, which corresponds to the principal symbol of ∆ g possessing symmetries. See the works of Cassanas [3] and Gornet [9] for other closely related work.
Let (M, g) be a compact Riemannian manifold with boundary satisfying the same assumptions as in [10]. In fact, all that is necessary is that one may construct a parametrix in the same form as [10] for the Neumann wave propagator. We assume that a compact Lie group G has a symplectic group action on T * M and it accounts for all the symmetries of the Hamiltonian p(x, ξ) = |ξ| 2 g , which is the principal symbol of ∆ g . That is, p(g.z) = p(z) for each z ∈ T * M and g ∈ G, where g.z or gz denotes the group action. The assumption is merely that G accounts for all such symmetries.
The symmetry in the Hamiltonian implies that the fixed point manifolds C T introduced in section 2 are in fact invariant under the group action: g.C T ⊂ C T for each g ∈ G. The main assumption we make is that the connected components of C T are the G-orbit of a particular periodic bicharacteristic γ and it has no other symmetry; i.e. assuming C T is already connected, C T = Gγ := {g.x; g ∈ G and x ∈ Im(γ)}. We write [γ] to denote this set, which is also a certain equivalence class of γ as described in section 2. We note that this assumption captures generic situations since on general manifolds without symmetry, C T has dimension 1 (see [8, p. 61]), and an increase in dimension should only come from a group symmetry.
Next, we set L = I − dΦ T (m) for m ∈ C T . In the context of group symmetries, the map I − dΦ T : T m S * M/ ker(L) → T m S * M/ ker(L) is in general not an isomorphism anymore [8,  will generically be an isomorphism [2, appendix A.1], [4,5,3]. The above map is what we will now refer to as the Poincaré map, which we denote I −P [γ] as described in section 2. We also assume the clean intersection hypothesis described in [8,10] so that C T is a submanifold and ker(L) = T m C T . See remark 4.3 for a geometric description of the clean intersection hypothesis in spherical symmetry.
The new part of the trace formula will come from ker(L 2 )/ ker(L) whose presence is already seen in the calculation of the normal Hessian for the stationary phase analysis in section 3.3. Based on the arguments in [4,5], there is a geometric scalar quantity, which we denote by d [γ] , defined on C T , that is nonvanishing and associated with the action of I −dΦ T on ker(L 2 )/ ker(L) [4, appendix]. This is a quantity directly related to periodic orbits in the reduced phase space formally written as T * M/G and defined in [1,Chapter 4]. Under certain geometric conditions, I − P [γ] and d [γ] will stay constant over C T = [γ]. For our SO(3) action, ker(L 2 )/ ker(L) was only one-dimensional determined by the infinitesimal generator of dilations in the dual variables. Geometrically, this corresponds to merely increasing the speed of a particular geodesic. Analytically, for a geodesic γ ν starting at ν ∈ S * M , it corresponds to the Jacobi field J(t) = tγ ν (t) (see [8, p. 70] for more details). Hence, in our SO(3) case, T m T * M/ ker(L 2 ) is indeed isomorphic to T m S * M/T m C T , but generally, ker(L 2 )/ ker(L) will be much bigger when there is symmetry.
Next, our assumptions imply that there is a natural surjective map [0, T γ ) × G → C T which usually fails to be injective, where T γ is the primitive period of γ. For example, in the analysis in section 3, one does not need the entire rotation group SO(3) to form C T , which would in fact create multiple copies, since certain elements of SO(3) already map γ to itself. If the map is locally injective, then following [5, Section 3], we assume there is a discrete subgroup I(γ) of G that carries γ to itself, and we denote N γ to be the number of elements in I(γ), corresponding to the analogous quantity appearing in proposition 2.2.
In summary, we have accounted for all the pieces in our trace formula (proposition 2.2) in a general Lie group framework so that our formula may be compared to the formula appearing in [5] albeit in a different setting. denote σ = ρ −1 . We will later choose ρ so that ρ(1) = 1 and ρ(R) > 0, which makes G ρ : A →B(0, 1) \ B(0, ρ(R)) a diffeomorphism.
For convenience, we change variable from r to s = σ(r). Since the integrand is strictly positive and C k , the function ρ is a strictly increasing C k+1 function as claimed. We also claimed earlier that ρ satisfies ρ(1) = 1 and ρ(R) > 0, and these properties can be read in the representation (C.4).
Remark C.2. The diffeomorphism of proposition C.1 is in fact radial if n ≥ 3. It will also be necessarily radial in two dimensions if the metric is invariant under the action of SO (2), but also that of O(2).
Example D.2. Consider the closed hemisphere of S 3 or any other S n , n ≥ 2. Now all broken rays are periodic, and all primitive periods are 2π. The primitive length spectrum is degenerated into a single point. The Herglotz condition fails at the boundary (but only there), and the endpoints of all maximal geodesics are conjugate. The fixed point set has full dimension; for n = 3 we have dim(C 2π ) = 5. This manifold is non-trapping.
Example D.3. Consider the closed hemisphere H ⊂ S 2 . In polar coordinates we can identify H near the boundary (equator) with (R, 1] × [0, 2π], where the angles 0 and 2π are identified in the obvious way. For any L ∈ (0, 2π), we can take the submanifold (R, 1]×[0, L] and identify angles 0 and L. If L / ∈ πQ, then no broken ray near the boundary is periodic. By angular rescaling we may write this as a rotation invariant metric near the boundary of the closed discD and continue it to a rotation invariant metric in the whole disc. Now there is an open set of points on the sphere bundle containing no periodic orbits, and this open set can be made large. Only rays going close enough to the center of the disc will contribute to the length spectrum. This manifold is also non-trapping, the Herglotz condition fails at the boundary, and the countable conjugacy condition is violated.