Spherically symmetric terrestrial planets with discontinuities are spectrally rigid

TL;DR

This paper proves spectral rigidity for spherically symmetric planets with interior discontinuities, allowing for different metrics and wave types, which helps in understanding planetary interior structures.

Contribution

It introduces a new spectral rigidity result for manifolds with interior metric discontinuities, applicable to planetary models with fluid cores and phase transitions.

Findings

Spectral rigidity holds for spherically symmetric manifolds with interior interfaces.

The length spectrum of the Euclidean ball is shown to be simple.

The proof utilizes a novel trace formula for these manifolds.

Abstract

We establish spectral rigidity for spherically symmetric manifolds with boundary and interior interfaces determined by discontinuities in the metric under certain conditions. Rather than a single metric, we allow two distinct metrics in between the interfaces enabling the consideration of two wave types, like P- and S-polarized waves in isotropic elastic solids. Terrestrial planets in our solar system are approximately spherically symmetric and support toroidal and spheroidal modes. Discontinuities typically correspond with phase transitions in their interiors. Our rigidity result applies to such planets as we ensure that our conditions are satisfied in generally accepted models in the presence of a fluid outer core. The proof is based on a novel trace formula. We also prove that the length spectrum of the Euclidean ball is simple.

Tables

	(A1)	(A2)	(A3)	(A4)
Trace formula	X	X		X
Length spectral rigidity			X	X
Spectral rigidity	X	X	X	X

Equations

$$\gamma (t_i)\in \Gamma ,(d{\iota }_{\Gamma }{)}^{\ast }(\gamma (t_i),\dot{\gamma }(t_i^{-}))=(d{\iota }_{\Gamma }{)}^{\ast }(\gamma (t_i),\dot{\gamma }(t_i^{+})),$$

$$t\mapsto (r(t),\theta (t))$$

$$L(r_0,p):={\int }_{r_0}^1\frac{1}{c(r^{\prime }{)}^2\beta (r^{\prime };p)}d{r}^{\prime }$$

$$\alpha (p)$$

$$T(p)=2N(p)L(p)=2N(p){\int }_{R^{\ast }}^{r_{k-1}}\frac{1}{c(r^{\prime }{)}^2\beta (r^{\prime };p)}d{r}^{\prime }$$

$$\alpha (p)=2N(p){\int }_{R^{\ast }}^{r_{k-1}}\frac{p}{(r^{\prime }{)}^2\beta (r^{\prime };p)}d{r}^{\prime }.$$

$$\underset{i\to \infty }{\lim }{\gamma }_i=\gamma$$

$$\kappa =a(\theta -{\theta }_0{)}^{1/2}+\mathcal{O}(\theta -{\theta }_0)$$

$$\phi (\theta )=b(\theta -{\theta }_0{)}^{1/2}+\mathcal{O}(\theta -{\theta }_0)$$

$$\alpha (\theta )-{\alpha }_0=c(\theta -{\theta }_0)+\mathcal{O}((\theta -{\theta }_0{)}^2)$$

$${\alpha }_0{b}^{-1}(\theta -{\theta }_0{)}^{-1/2}+\mathcal{O}((\theta -{\theta }_0{)}^{1/2})\in \mathbb{N}.$$

$$c(t-T(p)+i0{)}^{-k}{i}^{N(p)}n(p)Q(p)L(p){\left|p^{-2}{\partial }_p\alpha \right|}^{-1/2}$$

$$n(p_1)Q(p_1){\left|p_1^{-2}{\partial }_p\alpha (p_1)\right|}^{-1/2}\ne n(p_2)Q(p_2){\left|p_2^{-2}{\partial }_p\alpha (p_2)\right|}^{-1/2}.$$

$${\ell }^{\prime }(0)=c_0(r_1(0){)}^{-1}{r}_1^{\prime }(0).$$

$$h(r)={\int }_r^{1}f(s){\left[1-{\left(\frac{rc(s)}{sc(r)}\right)}^2\right]}^{-1/2}\frac{\mathrm{d}s}{c(s)}$$

$${\partial }_{\tau }\ell (\tau ,r^{\tau }){\upharpoonright }_{\tau =0}=2nh(r^{\prime }),$$

$$u={}_n{\mathbf{D}}_l{Y}_l^m,$$

$$\mathbf{D}=U(r)(-k^{-1})[-\hat{\theta }(\sin \theta {)}^{-1}{\partial }_{\psi }+\hat{\psi }{\partial }_{\theta }],$$

$$Y^{m}_{l}(\theta,\psi)=\left\{\begin{array}[]{rcl}\sqrt{2}X^{\left\lvert m\right\rvert}_{l}(\theta)\cos(m\psi)&\text{if}&-l\leq m<0,\\ X^{0}_{l}(\theta)&\text{if}&m=0,\\ \sqrt{2}X^{m}_{l}(\theta)\sin(m\psi)&\text{if}&0<m\leq l,\end{array}\right.$$

$$X_l^m(\theta )=(-{)}^m\sqrt{\frac{2l+1}{4\pi }}\sqrt{\frac{(l-m)!}{(l+m)!}}{P}_l^m(\cos \theta ),$$

$$P_l^m(\cos (\theta ))=(-{)}^m\frac{1}{2^{l}l!}(\sin \theta {)}^m{\left(\frac{1}{\sin \theta }\frac{\mathrm{d}}{\mathrm{d}\theta }\right)}^{l+m}(\sin \theta {)}^{2l}.$$

$$[-r^{-2}{\partial }_r{r}^2\mu {\partial }_r+r^{-2}{\partial }_r\mu r-r^{-1}\mu {\partial }_r+r^{-2}(-1+k^2)\mu ]U-{\omega }^2\rho U=0,$$

$$T(U)=\mathcal{N}U,\mathcal{N}=\mu {\partial }_r-r^{-1}\mu$$

$$U_{+}{\upharpoonright }_{r=b}$$

$${\mu }_{+}{\partial }_r{U}_{+}{\upharpoonright }_{r=b}$$

$$\{Y_l^{-l},\dots ,Y_l^l\}.$$

$${}_n{I}_l=I_d:={\int }_R^1{\left|U_d(r)\right|}^2\rho (r)dr$$

$$G(x,x_0,t)=\mathrm{Re}\sum _{l=0}^{\infty }\sum _{n=0}^{\infty }{}_n{\mathbf{D}}_l({}_n{\mathbf{D}}_l{)}_0\sum _{m=-l}^l{Y}_l^m(\theta ,\psi )Y_l^m({\theta }_0,{\psi }_0)\frac{e^{i{}_n{\omega }_{l}t}}{i({}_n{\omega }_l)({}_n{I}_l)}.$$

$$\sum _{m=-l}^l\iint Y_l^m(\theta ,\psi {)}^2\sin \theta d\theta d\psi =2l+1$$

$$\sum _{m=-l}^l\iint (-k^{-2}){\left|[-\hat{\theta }(\sin \theta {)}^{-1}{\partial }_{\psi }+\hat{\psi }{\partial }_{\theta }]Y_l^m(\theta ,\psi )\right|}^2\sin \theta d\theta d\psi =2l+1.$$

$${\int }_{M}G(x,x,t)\rho (x)dx=\sum _{l=0}^{\infty }\sum _{n=0}^{\infty }(2l+1)\mathrm{Re}\left\{\frac{e^{i{}_n{\omega }_{l}t}}{i({}_n{\omega }_l)}\right\}$$

Full text

Acknowledgements

MVdH was supported by the Simons Foundation under the MATH + X program, the National Science Foundation under grant DMS-1815143, and the corporate members of the Geo-Mathematical Imaging Group at Rice University. JI was supported by the Academy of Finland (projects 332890 and 336254). We thank Chunquan Yu for help with composing figure 2. We also thank Gabriel Paternain and Anthony Várilly-Alvarado for some helpful discussions. The authors greatly appreciate the detailed suggestions made by an anonymous referee which improved this paper.

1. Introduction

We establish spectral rigidity for spherically symmetric manifolds with boundary and interfaces determined by discontinuities in the metric. We study the recovery of a (radially symmetric Riemannian) metric or wave speed containing jump discontinuities along finitely many $C^{\infty}$ hypersurfaces. To our knowledge, it is the first such result pertaining to a manifold with boundary and a piecewise continuous metric.

Terrestrial planets in our solar system are approximately spherically symmetric. On the one hand, the deviation from such a symmetry becomes apparent only at high eigenfrequencies. On the other hand, our results provide a stable approximation upon truncating the spectrum of eigenfrequencies. Discontinuities arise largely due to phase transitions. Hence, their radial depths play an important role in determining the thermal structure and chemical composition of planets as well as the dynamics of their interiors [32]. The question of spectral rigidity is behind the validity of PREM [18] which is still widely used as a reference in linearized tomography. More interestingly, in space exploration such as the current NASA’s InSight mission to Mars [27], with a single data point, spectral data could provide the leading information about its interior; other missions are being proposed.

The results presented, here, are an extension of our previous result [14] where we proved a spectral rigidity for a smooth metric on a radial manifold. Allowing for certain discontinuities in the metric adds a new level of challenge for several reasons. First, the geodesics in such a manifold get reflected and transmitted when they hit an interface, creating a complex geometry for the analysis. In addition, we allow such geodesics to hit an interface at certain critical angles where a scattered ray can intersect an interface tangentially or “glide” along an interface. We also recover the location of the interfaces and do not assume that they are known.

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 ball $M=\bar{B}(0,1)\subset\mathbb{R}^{3}$, with the metric $g(x)=c^{-2}(\left\lvert x\right\rvert)e(x)$, where $e$ is the standard Euclidean metric and $c\colon(0,1]_{r}\to(0,\infty)$ is a function satisfying suitable conditions, where $r=|x|$ is the radial coordinate. We work in dimension three but our result on length spectral rigidity (Theorem 1.2) carries over to higher dimensions, and our methods to prove spectral rigidity (Theorem 1.6) may be generalized to higher dimensions.

We assume $c(r)$ has a jump discontinuity at a finite set of values $r=r_{1},\dots,r_{K}$; that is $\lim_{r\to r_{i}^{-}}c(r)\neq\lim_{r\to r_{i}^{+}}c(r)$ for each $i$. Our assumption is the smooth Herglotz condition: $\frac{\mathrm{d}}{\mathrm{d}r}(r/c(r))>0$ is satisfied everywhere away from the discontinuities of $c$, but we note that $c$ is allowed to either increase or decrease across an interface. We note that the natural extension of the Herglotz condition when $c$ is smooth to our case when $c$ has discontinuities is to view $c$ as a distribution and require $\frac{\mathrm{d}}{\mathrm{d}r}(r/c(r))>0$ in the distributional sense. If $c$ has a jump discontinuity at $r=r_{i}$, this distributional condition implies $\lim_{r\to r_{i}^{-}}c(r)>\lim_{r\to r_{i}^{+}}c(r)$. This would be too restrictive since radial models of Earth (PREM) and Mars (T13) (see [16]) satisfy the smooth Herglotz condition but not this stronger distributional Herglotz condition, since the jump across the core-mantle boundary differs in sign to the jumps at other interfaces. Hence, our smooth Herglotz condition is weaker to allow the jump across interfaces to have any sign. We also allow trapped rays that never interact with the boundary. Such rays just correspond to small but nonzero boundary amplitudes of modes. The assumption $\frac{\mathrm{d}}{\mathrm{d}r}(r/c(r))>0$ when $c$ is smooth is the Herglotz condition first discovered by Herglotz [23] and used by Wiechert and Zoeppritz [39].

By a maximal geodesic we mean a unit speed geodesic on the Riemannian manifold $(M,g)$ with each endpoint at the boundary $\partial M$ or an interface. 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$, and Snell’s law for geometric optics at the interfaces. If the initial and final points of a broken ray coincide at the boundary or an interface, 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 2.4) what will be called the countable conjugacy condition which ensures that up to rotation only countably many maximal geodesics have conjugate endpoints.

The length spectrum of a manifold $M$ with boundary is the set of lengths of all periodic broken rays on $M$. If the radial sound speed is $c$, we denote the length spectrum by $\operatorname{lsp}(c)$. We will introduce in Definition 2.3 the notion of closed basic rays, which are certain periodic rays that stay completely within a single layer. The set of lengths of such rays form the basic length spectrum $\operatorname{blsp}(c)$. We note that every broken ray is contained in a unique two-dimensional plane in $\mathbb{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, $\Delta_{c}=c^{3}\nabla\cdot c^{-1}\nabla$, on $M$ by $\operatorname{spec}(c)$, where we impose Neumann-type boundary conditions on both the inner and outer boundary. The spectrum $\operatorname{spec}(c)$ includes multiplicity, not just the set spectrum.

Some earlier results in tensor tomography, the methods of which are related to ours, may be found in [30, 31, 33, 34]. Let us now enumerate the various geometric assumption we make in this manuscript for easy reference.

1.1. Herglotz and other conditions

There are several geometric assumptions we make that we shall enumerate here:

1. (A1)

“Periodic conjugacy condition.” This is an analog of the clean intersection hypothesis used in [21, 17, 14]; (see Definition 2.5). 2. (A2)

“Principal amplitude injectivity condition.” This is an analog to assuming simplicity of the length spectrum. (see section 2.3). 3. (A3)

“Countable conjugacy condition” (Definition 2.4). 4. (A4)

Smooth Herglotz condition: $\frac{d}{dr}\frac{r}{c(r)}>0$ away from the discontinuities.

These assumptions allow us to prove that the singular support of the wave trace includes the basic length spectrum. Assumption (A1) is a standard assumption (normally referred to as the clean intersection hypothesis when $c$ is smooth) when calculating the trace singularity by a stationary phase method to ensure that the critical manifolds are non-degenerate and the phase function is Bott-Morse nondegenerate (see [17, 21]). A ubiquitous issue in computing a trace formula is the possibility of cancellations between the contributions of two components of the same length that are not time reversals of each other to the wave trace. One usually assumes “simplicity” of the length spectrum so that any two rays with a given period are either rotations of each other or time reversals of each other, but since our trace formula computation is more explicit, we have a slightly weaker assumption (A2) to take care of this issue. Assumptions (A1), (A2), and (A4) are needed for the trace formula (Proposition 4.1), and all four assumptions are needed for spectral rigidity (Theorem 1.6), while only assumptions (A3) and (A4) are used to prove length spectral rigidity (Theorem 1.2). Below, we provide a chart for easy reference regarding which assumptions are needed for each theorem:

*Genericity of the assumptions: * A key novelty of the present paper is that the final proof is essentially a length spectral rigidity result. On a general Riemannian manifold the very existence of closed orbits is a complicated matter [3, 4, 20], and having a boundary makes the problem even harder. In some geometries, such as negatively curved closed manifolds, periodic orbits are known to exist and indeed spectral rigidity has been proven on such manifolds [12, 11]. In the case of spherical symmetry we only need a very mild condition to ensure that there are enough periodic orbits and they behave cleanly enough; this is called the countably conjugacy condition. Thus, we only require assumptions (A3) and (A4) to show length spectral rigidity. We expect these condition to be generic on the class of manifolds we study, but verifying it is not simple check, so we leave it as a conjecture that assumptions (A1), (A2), and (A3) are generic among metrics that satisfy the Herglotz condition.

In Appendix A, we prove that the class of metrics satisfying all of our assumptions is nonempty and a genericity result. Specifically, we construct an uncountable family of metrics that satisfy assumptions (A.1), (A.3), and (A.4). For assumption (A.2), which is an analog of having a simple length spectrum in our setting and is more subtle to verify, we prove that within this class, assumption (A.3) is generically satisfied.

In Appendix A.2, we prove that the Euclidean disk does have a simple length spectrum, thereby satisfying all of the assumptions and demonstrating that the set of metrics satisfying assumptions (A.1)-(A.4) is nonempty. The proof relies on advanced concepts from number theory, even in this geometrically simplest case. This suggests that finding examples in non-Euclidean geometry will be difficult.

1.2. Main results

Here we present our main theorems, which follow a discussion of the notation we use for the geometry. Let $A(r^{\prime},r^{\prime\prime})=\bar{B}(0,r^{\prime\prime})\setminus B(0,r^{\prime})\subset\mathbb{R}^{3}$ be the closed annulus in a Euclidean space where $r^{\prime\prime}>r^{\prime}$. Fix $K\in\mathbb{N}$ and let $r_{k}\in(0,1)$ so that $1=:r_{0}>r_{1}>\cdots>r_{K}$. Assume $c(r)$ has jump discontinuities at each $r_{k}\in(0,1)$. Let $\Gamma=\bigcup_{k}\{r=r_{k}\}$ be the collection of interfaces together with $\partial M$, and denote $\Gamma_{k}:=\{r=r_{k}\}$. We sometimes refer to the smooth annular regions $A(r_{k},r_{k-1})$ as layers. We view $M$ as a Riemannian manifold with (rough) metric $g=c^{-2}dx^{2}$.

Definition 1.1**.**

Fix any $\varepsilon>0$ and $K\in\mathbb{N}$. We say that a collection of functions $c_{\tau}\colon[0,1]\to(0,\infty)$ indexed by $\tau\in(-\varepsilon,\varepsilon)$ is an admissible family of profiles if the following hold:

• •

There are radii $r_{k}\in(0,1)$ that depend $C^{1}$-smoothly on $\tau\in(-\varepsilon,\varepsilon)$ so that $1\eqqcolon r_{0}(\tau)>r_{1}(\tau)>\dots>r_{K}(\tau)>0$ for all $\tau\in(-\varepsilon,\varepsilon)$.

• •

For every $\tau\in(-\varepsilon,\varepsilon)$ the function $c_{\tau}$ is piecewise $C^{1,1}$ and satisfies the smooth Herglotz condition.

• •

The only singular points of each function $c_{\tau}$ are the radii $r_{k}(\tau)$ where it has a jump discontinuity.

• •

Within each annulus $A(r_{k}(\tau),r_{k-1}(\tau))$ the profile $c_{\tau}$ satisfies the countably conjugacy condition for all $\tau\in(-\varepsilon,\varepsilon)$.

• •

We assume that $(r,\tau)\mapsto c_{\tau}(r)$ is $C^{1}$ at all points where $r\notin\{r_{1}(\tau),\dots,r_{K}(\tau)\}$.

Recall from the introduction that the length spectrum of a manifold $M$ with boundary is the set of lengths of all periodic broken rays on $M$ and we denote the length spectrum by $\operatorname{lsp}(c)$. We will introduce in Definition 2.3 the notion of closed basic rays, which are certain periodic rays that stay completely within a single layer. The set of lengths of such rays form the basic length spectrum $\operatorname{blsp}(c)$.

Our main theorem provides the rigidity of the basic length spectrum in the presence of “countable noise”. Choosing the “noise” suitably gives corollaries for the full length spectrum. Missing or spurious points in the length spectrum or some amount of degeneracy do not matter. The “noise” can be of the same size as the data, and this will play a role in the case of multiple wave speeds.

Theorem 1.2**.**

Fix any $\varepsilon>0$ and $K\in\mathbb{N}$, and let $c_{\tau}(r)$ be an admissible family of profiles with discontinuities at $r_{k}(\tau)$ for all $k=1,\dots,K$.

Let $\operatorname{blsp}(\tau)$ denote the basic length spectrum of the ball $\bar{B}(0,1)$ with the velocity profile $c_{\tau}$. Suppose $\operatorname{blsp}(\tau)$ is countable for all $\tau$. Let $S(\tau)$ be any collection of countable subsets of $\mathbb{R}$ indexed by $\tau$. If $\operatorname{blsp}(\tau)\cup S(\tau)=\operatorname{blsp}(0)\cup S(0)$ for all $\tau\in(-\varepsilon,\varepsilon)$, then $c_{\tau}=c_{0}$ and $r_{k}(\tau)=r_{k}(0)$ for all $\tau\in(-\varepsilon,\varepsilon)$ and $k=1,\dots,K$.

The theorem has two immediate corollaries. The first one concerns the whole length spectrum, and the second one the length spectrum of two velocity profiles.

Corollary 1.3** (Length spectral rigidity of a layered planet with moving interfaces).**

Fix any $\varepsilon>0$ and $K\in\mathbb{N}$, and let $c_{\tau}(r)$ be an admissible family of profiles with discontinuities at $r_{k}(\tau)$ for all $k=1,\dots,K$. Suppose that the length spectrum for each $c_{\tau}$ is countable in the ball $\bar{B}(0,1)$.

Let $\operatorname{lsp}(\tau)$ and $\operatorname{blsp}(\tau)$ denote the length spectrum and the basic length spectrum of the ball $\bar{B}(0,1)$ with the velocity profile $c_{\tau}$. Suppose either one of the following holds:

• •

$\operatorname{lsp}(\tau)=\operatorname{lsp}(0)$* for all $\tau\in(-\varepsilon,\varepsilon)$.*

• •

$\operatorname{blsp}(\tau)=\operatorname{blsp}(0)$* for all $\tau\in(-\varepsilon,\varepsilon)$.*

Then $c_{\tau}=c_{0}$ and $r_{k}(\tau)=r_{k}(0)$ for all $\tau\in(-\varepsilon,\varepsilon)$ and $k=1,\dots,K$.

Corollary 1.4** (Length spectral rigidity with two polarizations).**

Fix any $\varepsilon>0$ and $K\in\mathbb{N}$, and let $c^{i}_{\tau}(r)$ with both $i=1,2$ be an admissible family of profiles with discontinuities at $r_{k}(\tau)$ for all $k=1,\dots,K$.

Consider all periodic rays which are geodesics within each layer and satisfy the usual reflection or transmission conditions at interfaces, but which can change between the velocity profiles $c^{1}_{\tau}$ and $c^{2}_{\tau}$ at any reflection and transmission. Suppose that the length spectrum of this whole family of geodesics, denoted by $\operatorname{lsp}(\tau)$, is countable in the ball $\bar{B}(0,1)$.

If $\operatorname{lsp}(\tau)=\operatorname{lsp}(0)$ for all $\tau\in(-\varepsilon,\varepsilon)$, then $c^{i}_{\tau}=c^{i}_{0}$ for both $i=1,2$ and $r_{k}(\tau)=r_{k}(0)$ for all $\tau\in(-\varepsilon,\varepsilon)$ and $k=1,\dots,K$.

The “noise” set $S(\tau)$ of Theorem 1.2 plays an important role. One metric is recovered at a time, and all rays that have one leg following the other metric or different legs in different layers are treated as noise.

The proofs of the corollaries are immediate:

• •

For Corollary 1.3, simply let $S(\tau)=\operatorname{lsp}(\tau)$.

• •

For Corollary 1.4, study the basic length spectra of the profiles $c^{1}(\tau)$ and $c^{2}(\tau)$ independently of each other and let again $S(\tau)=\operatorname{lsp}(\tau)$.

Remark 1.5**.**

Some variations of Theorem 1.2 and its corollaries hold true. One can introduce an impermeable core and work with a finite number of layers that do not exhaust the ball. One can choose to include or exclude rays with reflections from the lower boundary $r_{K}(\tau)$ and the results remain true for this smaller length spectrum, at least when $r_{K}$ is independent of $\tau$. The proofs are immediate adaptations of the one we give.

Recall the Neumann spectrum of the Laplace–Beltrami operator is denoted $\operatorname{spec}(c)$, where we impose Neumann-type boundary conditions (we can allow for other boundary conditions cf. section 4.2).

Theorem 1.6** (Spectral rigidity with moving interfaces).**

Fix any $\varepsilon>0$ and $K\in\mathbb{N}$, and let $c_{\tau}(r)$ be an admissible family of profiles with discontinuities at $r_{k}(\tau)$ for all $k=1,\dots,K$. Suppose that the length spectrum for each $c_{\tau}$ is countable in the ball $\bar{B}(0,1)\subset\mathbb{R}^{3}$. Assume also that the length spectrum satisfies the principal amplitude injectivity condition and the periodic conjugacy condition.

Suppose $\operatorname{spec}(\tau)=\operatorname{spec}(0)$ for all $\tau\in(-\varepsilon,\varepsilon)$. Then $c_{\tau}=c_{0}$ and $r_{k}(\tau)=r_{k}(0)$ for all $\tau\in(-\varepsilon,\varepsilon)$ and $k=1,\dots,K$.

Proof.

The spectrum determines the trace of the Green’s function by Proposition 4.1. As $\operatorname{spec}(\tau)=\operatorname{spec}(0)$ for all $\tau$, the trace is independent of $\tau$ and so are its singularities. The singularities are contained in the set $\operatorname{lsp}(\tau)$ by Proposition 4.1. We apply Theorem 1.2 to pass from length spectral information to geometric information.

We set $S(\tau)$ to be the singular support of the trace. Every length of a basic periodic broken ray only appears once in the whole length spectrum by assumption, whence there is a singularity for every basic length. Therefore $\operatorname{blsp}(\tau)\subset S(\tau)$. Now Theorem 1.2 implies the claim. ∎

Planets are full balls, but Theorem 1.6 holds for an annulus as well. Cf. Remark 1.5.

Remark 1.7** (Implications for planets).**

The theorem is stated for a scalar operator (the Laplace-Beltrami operator), but the proof extends to the radial elastic case and thus, round planets by considering the toroidal modes associated with the shear wave speed and their corresponding eigenfrequencies. The proof of the theorem is using a trace formula to recover the basic length spectrum from the spectrum and employ the length spectral rigidity results. See sections 4.1 and 4.2, where we initially start the proof of the trace formula using toroidal modes and show why the argument is identical for the scalar Laplace-Beltrami operator. In that case, we work inside an annulus with an inner boundary representing the core-mantle boundary for more generality. By considering toroidal modes, the argument for proving a trace formula for spheroidal modes that involve two wave speeds becomes more transparent and is discussed in section 4.4. Hence, by considering the spectrum of the radial isotropic elastic operator with natural boundary conditions, our arguments may be generalized to recover both elastic wave speeds using Corollary 1.4.

Remark 1.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. For spectral rigidity, we assume dimension three to ease the computation of the trace formula since it allows us to compute the leading order asymptotics of the eigenfunctions explicitly.

This paper will be essentially divided into parts. The first part is proving length spectral rigidity. In the second part, we prove the trace formula in our setting, and as a corollary, we prove the spectral rigidity theorem.

1.3. Reasonableness of radial models

Spherically symmetric Earth models are widely used in geophysics and there are a number of results showing how well such models fit seismic data. The $P$ and $S$ wave speeds are denoted $c_{P}$ and $c_{S}$. There are several important questions to address when using PREM to analyze seismic data.

Question 1. What is the uncertainty in the best-fitting spherical average profile?

The classic reference for this question is Lee and Johnson in [26]. They suggest that the extremal bounds in the upper mantle are around 0.6 km/s (around 6 %) for $c_{P}$ and 0.4 km/s for $c_{S}$ (around 7 %). In the lower mantle, it is around 0.18 km/s (around 2 %) for $c_{P}$, and 0.14 km/s (around 2 %) for $c_{S}$. Note that the bounds increase in the lowermost mantle and especially in the crust.

Question 2. What is the standard deviation of the residuals to the spherical average model, as a function of depth?

In theory, residuals can be calculated as a function of depth for any global tomographic model. However, this information is not always presented. A good, thorough, recent example is the SP12RTS model [25]. Their figure 9a shows that variations are smallest in the mid-mantle (standard deviations of around 0.1 % for $c_{P}$, 0.2 % for $c_{S}$) and increase towards the surface (to around 1.0 % for both $c_{P}$ and $c_{S}$) and towards the CMB (to around 0.3 % for $c_{P}$, and 0.5 % for $c_{S}$).

Question 3. What is the measurement uncertainty in the wave speed at a given point in a typical tomographic model?

Very few groups have given robust estimates of point-wise measurement uncertainties, and the best study to date could be the Bayesian study by Burdick and Lekić in [7]. They find the standard deviation in estimates of 0.25 % $dc_{P}/c_{P}$ (so, for example the anomaly in California at 10 km depth might be 1.00 % +/- 0.25 %). We are not aware of any similar estimates for $c_{S}$, but they would most likely be more uncertain.

Question 4. In a given region, what is the typical variation in the absolute wavespeed?

Near Earth’s surface, there are huge lateral variations in wavespeed, for example between continental and oceanic regions (for example, at a depth of 50 km, mountain belt may have a $c_{P}$ of 6.1 km/s, while an ocean basin may have a $c_{P}$ of 8.1 km/s at the same radial coordinate, a variation of 25 %. However, within a given region type (e.g. ’island arc’ or ’mountain belt’), typical variations around 0.3 km/s for $c_{P}$ (an authoritative reference is [28]; see their fig. 3b), which is about 5 %. Variations in $c_{S}$ can be larger because $c_{S}$ is more strongly affected by fluids and temperature (partial melting and anelasticity). The reference given does not address $c_{S}$.

2. Unraveling assumptions

Let us give the relevant definition and assumptions on the geometry of the problem. Recalling from the previous section, fix $K\in\mathbb{N}$ and let $r_{k}\in(0,1)$ so that $1=:r_{0}>r_{1}>\cdots>r_{K}$. Assume $c(r)$ has jump discontinuities at each $r_{k}\in(0,1)$. Let $\Gamma=\bigcup_{k}\{r=r_{k}\}$ be the collection of interfaces together with $\partial M$, and denote $\Gamma_{k}:=\{r=r_{k}\}$. We view $M$ as a Riemannian manifold with (rough) metric $g=c^{-2}dx^{2}$. We showed in [14] that any rotation symmetric Riemannian manifold with the Herglotz condition is of this form. The same is true in the presence of jumps with essentially the same proof we used in the smooth setting.

2.1. Geodesics in a spherically symmetric model with interfaces

On the three-dimensional manifold $M$ the phase space of the unit speed geodesic flow has dimension $5$. Due to rotation symmetry most of these dimensions are superfluous, and the dimension of the reduced phase space needed to represent all geodesics up to isometries of the manifold is only $2$. The dimension of the “reduced phase space” is $2$ for any ambient dimension $2$ or higher.

Two natural coordinates in this space are the radius $r$ (Euclidean distance to the origin) and the angular momentum denoted as $p$. Any geodesic is either radial or is contained in a unique plane through the origin, so it suffices to study geodesics in $2$-dimensional disks. In dimension two, points on the disk can be described with polar coordinates $(r,\theta)$, and a geodesic $\gamma$ can be parameterized as $t\mapsto(r(t),\theta(t))$. We then have the explicit formula $p=p_{\gamma}=c(r(t))^{-2}r(t)^{2}\theta^{\prime}(t)$. The angular momentum (often called the ray parameter associated to $\gamma$) $p$ is conserved, even across discontinuities in the metric. Therefore trajectories of the geodesic flow in the $(r,p)$-plane are horizontal lines.

Much of the geometry is conveniently encoded in the function $\rho(r)=r/c(r)$. At a turning point (where $\dot{r}=0$) we have $\left\lvert p\right\rvert=\rho(r)$, and elsewhere $\left\lvert p\right\rvert<\rho(r)$. Therefore the reduced phase space is the subgraph of the function $\rho\colon(0,1]\to(0,\infty)$. The classical Herglotz condition states that $\rho^{\prime}(r)>0$ for all $r$. Three examples are given in figure 1.

Definition 2.1**.**

A (unit-speed) broken geodesic or ray in $(M,g)$ is a continuous, piecewise smooth path $\gamma:\mathbb{R}\supset I\to M$ such that each smooth piece is a unit-speed geodesic with respect to $g_{c}$ on $M\setminus\Gamma$, intersecting the interfaces $\Gamma$ at a discrete set of times $t_{i}\in I$. Furthermore, at each $t_{i}$, if the intersection is transversal, then Snell’s law for reflections and refraction of waves is satisfied. More precisely, a broken geodesic (parameterized by a time variable) can be written as $\gamma:(t_{0},t_{1})\cup(t_{1},t_{2})\cup\cdots\cup(t_{k-1},t_{k})\to M\setminus\Gamma$, which is a sequence of geodesics concatenated by reflections and refractions obeying Snell’s law: for $i=1,\dots,k-1,$

[TABLE]

where $\iota_{\Gamma}:\Gamma\to M$ is the inclusion map and $\dot{\gamma}(t_{i}^{\mp})=\lim_{t\to t_{i}^{\mp}}\gamma(t)$. Each restriction $\gamma\restriction_{(t_{i},t_{i+1})}$ is a maximal smooth geodesic that we call a leg of $\gamma$. For each $i$, note that $\gamma(t_{i})\in\Gamma_{k_{i}}$ for some $k_{i}$. One can view $\gamma$ as a concatenation of all of its legs. A leg $\gamma\restriction_{(t_{i},t_{i+1})}$ is reflected if the inner product of $\dot{\gamma}(t_{i}^{+})$ and $\dot{\gamma}(t_{i}^{-})$ with a normal vector to $\Gamma_{k_{i}}$ have opposite signs. If they have the same sign, it is a transmitted leg. If $\dot{\gamma}(t_{i}^{+})$ and $\dot{\gamma}(t_{i}^{-})$ are equal, then $\gamma\restriction_{(t_{i-1},t_{i+1})}$ is a grazing leg or ray; in this case, $\dot{\gamma}(t_{i}^{\pm})$ is tangent to $\Gamma$. The only other situation is when $\dot{\gamma}(t_{i}^{+})$ is tangent to $\gamma$ while $\dot{\gamma}(t_{i}^{-})$ is not (or vice versa); in this case $\gamma\restriction_{(t_{i},t_{i+1})}$ is called a gliding ray or leg because it travels along $\Gamma_{i_{k}}$. A ray with no gliding legs will be called a non-gliding ray. Our results will also extend to the elastic setting, which has two wave speeds $c_{P}$ and $c_{S}$ corresponding to pressure waves and shear waves. In this case, the definition of broken rays is identical except that each leg can either be a geodesic with the metric $g_{c_{P}}$ or $g_{c_{S}}$.

We follow the discussion and notation in [14, Section 2.1]. Assume for the moment $n=2$ since due to spherical symmetry, rays are confined to a disk, and equip the annulus $M=A(1,r)$ with spherical coordinates $\theta,r$. Fix a broken geodesic $\gamma$ whose endpoints are both located at a particular interface $r_{i}$ for some $i\in\{0,\dots,K\}$. We denote $\alpha=\alpha(p)$ to be the epicentral distance between both endpoints of $\gamma$, where $p=p_{\gamma}$ is the ray parameter associated to $\gamma$. It is the angular distance that $\gamma$ travels. It may happen that $\alpha(p)>2\pi$ if the geodesic winds around the origin several times.

Each leg can be parameterized as

[TABLE]

over some maximal interval $I$ associated to the leg. Using both of the conserved quantities $c(r(t))^{-2}[r^{\prime}(t)^{2}+r(t)^{2}\theta^{\prime}(t)^{2}=1$ and $p=c(r(t))^{-2}r(t)^{2}\theta^{\prime}(t)$ (angular momentum) we can compute $\alpha_{\gamma}$ explicitly following [14, Equation (2.2)].

Let $R^{*}$ be the smallest radius that $\gamma$ passes through, and there is a unique $k$ such that $r_{k}\leq R^{*}<r_{k-1}$. We refer to $R^{*}$ as the radius of $\gamma$ and it may coincide with an interface or a boundary. Next, $\gamma$ will have a certain number of legs in each of the annular regions $A(r_{k-1},r_{k})),A(r_{k-2},r_{k-1}),\dots A(r_{0},r_{1})$. Since $\gamma$ might stay just within a single (or more) annular region, there could be zero legs in one or more of the annuli. By definition of $R^{*}$, $\gamma$ has no legs in $A(r_{k},r_{K})$. We denote $n_{j}$ to be half of the number of legs of $\gamma$ in $A(r_{j-1},r_{j})$. Next we introduce a certain quantity $\beta^{2}:=c(r)^{-2}-r^{-2}p^{2}$.

Analogous to [14], the length of a broken geodesic with only transmitted legs, starting in $r=r_{0}$ and ending at $r=1$ is an integer multiple of the quantity

[TABLE]

If $r_{0}=R^{*}$ is the radius of $\gamma$, then $R^{*}$ is a function of $p$ and we will write $L(p)$. With this notation and using the computation for epicentral distance in [14], one can also find an explicit formula for $\alpha_{\gamma}(r):$

[TABLE]

Definition 2.2**.**

Following Hron in [24], those waves which travel from the source to the receiver along different paths but with identical travel-times are kinematically equivalent and are called kinematic analogs. We will refer to two different rays connecting source and receiver with the same ray parameter and travel time as kinematic analogs. The groups of kinematic analogs may be further divided into subgroups of waves whose amplitude curves are identical. The members of this subgroup of phases may be called dynamic analogs. A sufficient condition for kinematic equivalence of two different broken rays $\gamma_{1}$ and $\gamma_{2}$ is they must have an equal number of legs in each layer along their paths. Since $\alpha(p_{\gamma})$ just measures the epicentral distance between the endpoints, $\alpha(p_{\gamma})$ will be the same for $\gamma$ and all of its kinematic analogs. We will say two non-gliding rays connecting source and receiver are dynamic analogs if they have the same ray parameter, travel time, and inside each $A(r_{k},r_{k-1})$, they have the same number of legs that are reflections starting at $\Gamma_{k}$, transmission starting at $\Gamma_{k}$, reflections starting at $\Gamma_{k-1}$ and transmissions starting at $\Gamma_{k-1}$. This is a sufficient condition to ensure that the principal amplitudes of the corresponding waves are identical. See [24] for examples and figures of kinematic and dynamic analogs.

For length spectral rigidity, we only require what we term basic closed rays.

Definition 2.3** (Basic rays).**

A broken ray is called basic if either it stays within a single layer and all of its legs are reflections from a single interface (type 1), or it is a radial ray contained in a single layer (type 2). A radial ray is defined to be a ray with zero epicentral distance. It will necessarily reflect from two interface and cannot be type 1. The first type of basic rays are analogs to the turning rays in [14] that formed $\operatorname{lsp}(c)$ in the notation there. A closed basic ray of the first type will be periodic, stay within a singular layer, and only consists of reflected legs from a single interface. We have illustrated basic and other periodic rays in Figure 2.

The lengths of periodic basic will suffice to prove length spectral rigidity so we define $\operatorname{blsp}(c)$ as the set of lengths of all periodic basic rays.

Computing the length and epicentral distance of basic rays is much simpler. Let $\gamma$ be a basic ray with radius $R^{*}$, ray parameter $p$, and lies inside inside $A(r_{k-1},r_{k})$. Then there is a unique $N(p)\in\mathbb{N}$ so that the length, denoted $T(p)$, of $\gamma$ is

[TABLE]

and

[TABLE]

Definition 2.4**.**

Consider geodesics in an annulus $A(a,b)$ equipped with a $C^{1,1}$ sound speed $c\colon(a,b]\to(0,\infty)$. We say that $c$ satisfies the countable conjugacy condition if there are only countably many radii $r\in(a,b)$ so that the endpoints of the corresponding maximal geodesic $\gamma(r)$ are conjugate along that geodesic.

We will only need the countable conjugacy condition with each layer, so we do not need a definition in the presence of discontinuities. We point out that “countable” includes the possibility that the set be empty or finite. Definition 2.4 is the same as the one given in [14].

We need an analog to the clean intersection hypothesis used in [14, 21] to prove a trace formula that also makes sure that the phase function is Bott-Morse nondegenerate when applying a stationary phase argument.

Definition 2.5**.**

We say that the radial wave speed $c$ satisfies the periodic conjugacy condition if for each periodic, nongliding ray with a ray parameter $p$, $\partial_{p}\alpha(p)\neq 0$. This condition ensures that the phase function in the stationary phase argument for computing the trace formula is Bott-Morse nondegenerate.

2.2. Gliding rays as limits

Consider a periodic broken ray $\gamma_{0}$ with a gliding leg of positive length.

We assume that gliding occurs at only one interface; this is ensured by the smooth Herglotz condition. We may rearrange the legs of the periodic broken ray without changing its length or essential geometry so that there is only one gliding leg per period.

We will argue that there is a sequence of periodic non-gliding broken rays $\gamma_{i}$ so that $\gamma_{i}\to\gamma_{0}$. This is very simple for any finite segment of a gliding broken ray; the subtlety lies in ensuring periodicity of the approximating rays. We will prove the following lemma.

Lemma 2.6**.**

Let $\gamma_{0}$ be a periodic broken ray with a gliding leg of positive length as described above. Then there is a sequence $\{\gamma_{i}\}_{i=1}^{\infty}$ of periodic, non-gliding broken rays such that

[TABLE]

Proof.

Let $x$ and $y$ be the final and initial point, respectively, of the gliding leg of $\gamma_{0}$, and let $\theta_{0}$ be the angle between $\gamma_{0}$ and the interface. We wish to find angles $\theta_{i}>\theta_{0}$ with the correct approximating and periodicity properties.

For any angle $\theta>\theta_{0}$, let the angle between the interface and the leg of the refracted ray in the lower layer be denoted by $\kappa$. In the limiting case $\kappa_{0}=0$ as the ray $\gamma_{0}$ glides on the interface. It follows from Snell’s law and a calculation that

[TABLE]

for some constant $a>0$.

When $\theta$ is slightly above $\theta_{0}$ — or when $\kappa>0$ is small — the opening angle of a single short diving leg under the interface is denoted by $\varphi(\theta)$. A simple calculation shows that $\varphi(\theta)$ is asymptotically comparable to $\kappa$, whence

[TABLE]

for some constant $b>0$.

Let the angle between the points $y$ and $x$ be $\alpha_{0}>0$. Starting from the point $x$ and following the broken ray near $\gamma_{0}$ with the initial angle $\theta\approx\theta_{0}$ we get a map $\theta\mapsto y(\theta)$. This map is $C^{1}$. Denote the angle between $y(\theta)$ and $x$ by $\alpha(\theta)$. This map is well defined in a neighborhood of $\theta_{0}$, as the relevant broken ray stays above the interface and total internal reflection is not an issue.

If $\alpha^{\prime}(\theta_{0})=0$, then the points $x$ and $y$ are conjugate along the non-gliding part of the broken ray $\gamma_{0}$. But this turns out not to be an issue.

Denoting $\alpha^{\prime}(\theta_{0})=c$, we have

[TABLE]

due to simple Taylor approximation.

We want to choose the angle $\theta>\theta_{0}$ so that an integer amount of these short diving legs connect $y(\theta)$ to $x$. The condition is $\alpha(\theta)/\varphi(\theta)\in\mathbb{N}$. Combining with equations (2.4), (2.5), and (2.6), we end up with the condition that

[TABLE]

Here the error term depends continuously on $\theta$, so the left-hand side of equation (2.7) obtains integer values infinitely many times as $\theta\to\theta_{0}+$. This gives us a choice of directions $\theta_{i}$ starting at $x$, and thus a sequence of periodic broken rays $\gamma_{i}$ which converge to $\gamma_{0}$.

This concludes the argument that every periodic broken ray with a gliding leg can be approximated by periodic non-gliding rays. ∎

2.3. Principal amplitude injectivity condition

We also need an assumption similar to “simplicity” of the length spectrum modulo the group action in order to recover the length spectrum when there are multiple components in the length spectrum. For a closed ray $\gamma$, denote $[\gamma]$ to be the equivalence class to include all rotations and dynamic analogs of $\gamma$ along with its time reversal. We will see that $[\gamma]$ has a particular contribution to the trace formula.

The principal contribution of $[\gamma]$ with ray parameter $p$ to the trace formula has the form (see (4.20))

[TABLE]

where $c$ is independent of $\gamma$, $Q(p)$ is a product of reflection and transmission coefficients, and $T(p)$ is the length of $\gamma$. Theoretically, there may be another class $[\gamma^{\prime}]$ with an identical period whose principal contribution to the trace cancels with that of $[\gamma]$, thereby preventing recovery of $T$.

We say that the length spectrum satisfies the principal amplitude injectivity condition if given two closed rays $\gamma_{1}$ and $\gamma_{2}$ with the same period and disjoint equivalence classes (so they must have different ray parameters $p_{1}$ and $p_{2}$), then

[TABLE]

We assume that $\operatorname{lsp}(c)$ satisfies the principal amplitude injectivity condition in order to prove Theorem 1.6.

2.4. Spherical symmetry

In section 1.3 we saw that spherical symmetry is a good approximation for the Earth. This symmetry is of tremendous technical convenience. The geodesic flow is integrable with simple conserved quantities (an orbital plane and an angular momentum) and many of our calculations can be done explicitly.

The geometry of periodic broken rays is poorly understood outside symmetric situations. It is not clear whether there are densely many such rays on a general manifold with boundary, nor whether the periodic rays are stable under deformations of the geometry.

On general manifolds, small smooth perturbations of a smooth metric only have a second order effect on the direction of the geodesics. However, small smooth deformations of an interface have a first order effect, and this increased sensitivity substantially complicates matters. Radial deformations of radial models are better behaved in that the conserved symmetry and deformed conserved quantities make the deformations tractable.

3. Proofs: Length spectral rigidity

3.1. Auxiliary results

We denote by $A(r_{1},r_{0})=\bar{B}(0,r_{1})\setminus B(0,r_{0})\subset\mathbb{R}^{n}$ the closed annulus in a Euclidean space.

Lemma 3.1**.**

Fix any $\varepsilon>0$ and $r_{1}\in(0,1)$, and any finite set $F\subset(0,1)$. Let $r(\tau)\in(0,1)$ depend $C^{1}$-smoothly on $\tau$. Let $c_{\tau}$ with $\tau\in(-\varepsilon,\varepsilon)$ be $C^{1,1}$ functions $[r_{1},1]\to(0,\infty)$ satisfying the Herglotz condition and the countable conjugacy condition and depending $C^{1}$-smoothly on $\tau$.

If $\partial_{\tau}c_{\tau}(r)\restriction_{\tau=0}\neq 0$ for some $r\in(r_{1},1)$, then there is a periodic broken ray $\gamma_{\tau}$ with respect to $c_{\tau}$ so that

• •

$\tau\mapsto\ell_{\tau}(\gamma_{\tau})$* is $C^{1}$ on $(-\delta,\delta)$ for some $\delta\in(0,\varepsilon)$,*

• •

$\partial_{\tau}\ell(\gamma_{\tau})\restriction_{\tau=0}\neq 0$, and

• •

the depth (minimum of Euclidean distance to the origin) of $\gamma_{0}$ is not in $F$.

Here $\ell_{\tau}$ is the length functional corresponding to the velocity profile $c_{\tau}$.

While in our application we have $F=\emptyset$, we include this freedom in the lemma so that finitely many problematic depths can be avoided if needed.

We say that a broken ray is radial if it is contained in a one-dimensional linear (not affine) subspace of $\mathbb{R}^{n}$.

Lemma 3.2**.**

Fix any $\varepsilon>0$. Let $c_{\tau}\colon(0,1]\to(0,\infty)$ be a family if $C^{1,1}$ functions depending smoothly on $\tau\in(-\varepsilon,\varepsilon)$. Let $r(\tau)\colon(-\varepsilon,\varepsilon)\to(0,1)$ be $C^{1}$.

Let $\ell_{\tau}$ be the length of the radial geodesic between $r=r_{1}$ and $r=1$. If $\partial_{\tau}c_{\tau}(r)\restriction_{\tau=0}=0$ for all $\tau\in(-\varepsilon,\varepsilon)$, then

[TABLE]

3.2. Proof of Theorem 1.2

The idea of the proof is as follows: We first show that $c_{\tau}$ is independent of $\tau$ within the first layer. Then we show that the first interface is also independent of $\tau$. After these steps we can “peel off” the top layer and repeat the argument for the second one. Countability of the basic length spectrum provides sufficient decoupling between the layers and between the “data” $\operatorname{blsp}(\tau)$ and the “noise” $S(\tau)$.

We give most arguments at $\tau=0$ first for definiteness, but the exact value of the parameter is unimportant.

Proof of Theorem 1.2.

Let us denote $f_{\tau}(r)=\partial_{\tau}c_{\tau}(r)$ and $\hat{S}(\tau)=\operatorname{blsp}(\tau)\cup S(\tau)$.

Take any $r\in(r_{1}(0),1)$. If $\partial f_{0}(r)\neq 0$, then by Lemma 3.1 there is a family of basic periodic broken rays $\gamma_{\tau}$ for which the length map $\tau\mapsto\ell(\gamma_{\tau})$ is $C^{1}$ in a neighborhood of $\tau=0$ and the derivative at $\tau=0$ is non-zero.

As $\ell(\gamma_{\tau})\in\hat{S}(\tau)$ and by assumption $\hat{S}(\tau)=\hat{S}(0)$ for all $\tau$, this implies that the set $\hat{S}(0)$ contains a neighborhood of $\ell(\gamma_{0})$. This is in contradiction with countability of $\hat{S}(\tau)$, and so $\partial f_{0}(r)\neq 0$ is impossible.

We conclude that $\partial f_{\tau}(r)=0$ for all $r\in(r_{1}(0),1)$. The same argument can be repeated at any value of the parameter $\tau$, leading us to conclude that $\partial f_{\tau}(r)=0$ whenever $r\in(r_{1}(\tau),1)$.

If $r_{1}^{\prime}(0)\neq 0$, then by Lemma 3.2 the radial broken rays (which are basic and periodic with period twice their length) there is a family of broken rays whose lengths vary differentiably and with a non-zero derivative at $\tau=0$. This contradicts countability as above. The same argument is valid for any value of $\tau$, so we conclude that $r_{1}^{\prime}(\tau)=0$ for all $\tau\in(-\varepsilon,\varepsilon)$.

We have thus found that $r_{1}(\tau)=r_{1}(0)$ and $c_{\tau}(r)=c_{0}(r)$ for all $\tau\in(-\varepsilon,\varepsilon)$ and $r\in(r_{1}(0),1)$. We may now turn our attention to the annulus $A(r_{2}(\tau),r_{1}(\tau))$, whose top interface is now fixed at $r=r_{1}(0)=r_{1}(\tau)$ for all $\tau$. Repeating the same argument in this annulus shows that both the velocity profile in this annulus and the location of the second interface are independent of $\tau$. Carrying on inductively, we exhaust all layers of the ball and find that the claim does indeed hold true. ∎

3.3. Proofs of the lemmas

Lemma 3.1 is a small variation of the reasoning in [14], rewritten in a way that is useful in the presence of interfaces. The proof is concise; the reader is invited to refer to [14] for details.

Proof of Lemma 3.1.

Consider the velocity profile for any fixed $\tau$. A maximal broken ray without reflections from the inner boundary component is determined uniquely up to rotation by its deepest point. Let us denote the essentially unique geodesic of depth $r\in(0,1)$ by $\gamma_{r}^{\tau}$. For a subset $P^{\tau}\subset(r_{1},1)$ the corresponding broken rays are periodic, and we denote the minimal period by $\ell(\tau,r)$.

A periodic broken ray with respect to $c_{0}$ is called stable if there is $\delta\in(0,\varepsilon)$ so that there is a family of paths $\gamma^{\tau}\colon{\mathbb{R}}\to A(1,r_{1})$ which is $C^{1}$ in $\tau$ (and only continuously at reflection points) and each $\gamma^{\tau}$ is a periodic broken ray with respect to $c_{\tau}$. When such a family exists, let us denote the depth corresponding to the parameter $\tau\in(-\delta,\delta)$ by $r^{\tau}$. Let us denote by $C^{0}\subset P^{0}\subset(r_{1},1)$ the set of depths of stable periodic broken rays. It was shown in [14] that under the countable conjugacy condition and the Herglotz condition the set $C^{0}$ is dense in $[r_{1},1]$. Thus also $C^{0}\setminus F$ is dense.

Let us denote $f(r)=\partial_{\tau}c_{\tau}(r)\restriction_{\tau=0}$. Suppose that $f(r)\neq 0$ for some $r\in(r_{1},1)$. Due to the injectivity of generalized Abel transforms, the function

[TABLE]

is also non-trivial. As $h$ is continuous and $C_{0}$ is dense, there is $r^{\prime}\in C_{0}\setminus F$ so that $h(r^{\prime})\neq 0$.

The length $\ell(\tau,r^{\tau})$ of the family of periodic broken rays is differentiable in $\tau$ near $\tau=0$ because $r^{\prime}\in C^{0}$ and

[TABLE]

where $n$ is the (constant) winding number of the minimal period of $\gamma^{\tau}$.

Therefore the claimed derivative is indeed non-zero. ∎

The proof of Lemma 3.2 is a straightforward calculation and the statement is geometrically intuitive, so we skip the proof. The essential statement concerns simply the derivative of the length of a geodesic with respect to its endpoint.

4. The Trace formula and its proof

As in [14], we will prove a trace formula in order to recover part of the length spectrum, and then use the argument in the previous sections on length spectral rigidity in order to prove Theorem 1.6. Although the main theorems as stated in subsection 1.2 refer to the scalar operator $\Delta_{c}$, for greater generality, we initially consider the toroidal modes corresponding to the isotropic elastic operator (see [13, 14] for definitions). As in [14], the proof is identical when considering the scalar Laplace-Beltrami operator. This allows us to naturally consider and extend our results to spheroidal modes in section 4.4 where two waves speed are present. First, we give the general setup and state the trace formula as Proposition 4.1, followed by its proof.

4.1. Toroidal modes, eigenfrequencies, and trace formula

We now use spherical coordinates $(r,\theta,\psi)$. 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. Analytically, these eigenfunctions admit a separation in radial functions and real-valued spherical harmonics, that is,

[TABLE]

where

[TABLE]

in which $k=\sqrt{l(l+1)}$ 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}\mathbf{D}_{l}$ as the multiplication with ${}_{n}U_{l}(-k^{-1})$. In the above, $Y^{m}_{l}$ are spherical harmonics, defined by

[TABLE]

where

[TABLE]

in which

[TABLE]

The function, $U$ (a component of displacement), satisfies the equation

[TABLE]

where $\mu=\mu(r)$ is a Lamé parameter and $\rho=\rho(r)$ is the density, both of which are smooth, and $c=\sqrt{\mu/\rho}$. Also, $\omega={}_{n}\omega_{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

[TABLE]

which vanishes at the boundaries (Neumann condition). The transmission conditions are that $U$ and $T(U)$ remain continuous across the interfaces. If $r=b$ is an interface and $U_{\pm}$ represent two solutions on opposite sides of the interface, then in the high frequency limit as $\omega\to\infty$, the transmission conditions will amount to

[TABLE]

for the principal terms in the WKB expansion of the solution.

The radial equations do not depend on $m$ and, hence, every eigenfrequency is degenerate with an associated $(2l+1)$-dimensional eigenspace spanned by

[TABLE]

Following [40], let $d$ indicate the overtone number $n$ and the angular degree $l$. The radial eigenfunction $U_{d}(r)$ is independent of the order $m$. We define the inner product of the eigenfunctions:

[TABLE]

We use spherical coordinates $(r_{0},\theta_{0},\psi_{0})$ for the location, $x_{0}$, of a source, and introduce the shorthand notation $({}_{n}\mathbf{D}_{l})_{0}$ for the operator expressed in coordinates $(r_{0},\theta_{0},\psi_{0})$. We now write the (toroidal contributions to the) fundamental solution as a normal mode summation

[TABLE]

On the diagonal, $(r,\theta,\psi)=(r_{0},\theta_{0},\psi_{0})$ and, hence, $\Theta=0$. Here $\Theta$ is the angular epicentral distance, We observe the following reductions in the evaluation of the trace of (4.8):

• •

We will not normalize $U(r)$. Meanwhile, the spherical harmonic terms satisfy

[TABLE]

(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

[TABLE]

From the reductions above, we obtain

[TABLE]

or

[TABLE]

Let us also denote $\Sigma=\text{singsupp}(\operatorname{Tr}(\partial_{t}G))\subset\mathbb{R}_{t}.$

4.2. Connection between toroidal eigenfrequencies, spectrum of the Laplace–Beltrami operator, and the Schrödinger equation

We repeat the discussion in [14] to relate the spectrum of a scalar Laplacian, the eigenvalues associated to the vector valued toroidal modes, and the trace distribution $\sum_{l=0}^{\infty}\sum_{n=0}^{\infty}\ (2l+1)\cos(t{}_{n}\omega_{l})$.

We note that (4.3) and (4.4) for $U$ ensure that $v=UY^{m}_{l}$ satisfies

[TABLE]

where $P_{0}=r^{-1}(\partial_{r}\mu)$ is a [math]th order operator, $\omega^{2}$ is a particular eigenvalue, and $\mathcal{N}$ is as in (4.4). Hence $UY^{m}_{l}$ are scalar eigenfunctions for the self-adjoint (with respect to the measure $\rho\operatorname{d}\!x$) scalar operator $P$ with Neumann boundary conditions (on both boundaries) expressed in terms of $\mathcal{N}$.

The above argument shows that we may view the toroidal spectrum $\{{}_{n}\omega^{2}_{l}\}_{n,l}$ as also the collection of eigenvalues $\lambda$ for the boundary problem on scalar functions (4.13). Thus (4.12) can be written in the form

[TABLE]

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 $\operatorname{Tr}(\partial_{t}G)$, which is equal to $\sum_{\lambda\in\operatorname{spec}(P)}\cos(t\sqrt{\lambda})$ by the normalizations we have chosen.) Up to principal symbols, $P$ coincides with the $\Delta_{c}=c^{3}\nabla\cdot c^{-1}\nabla$ upon identifying $c^{2}$ with $\rho^{-1}\mu$. This means that the length spectra of $P$ and $\Delta_{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.

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,\theta)$ (which are spherical coordinates in dimension three by treating $\theta$ as coordinates on the $2$-sphere),

[TABLE]

where $\Delta_{\theta}$ is the Laplacian on the $2$-sphere.

Let us now simplify the PDE (4.13) for $v$. Let $Y(\theta)$ be an eigenfunction of $\Delta_{\theta}$ with eigenvalue $-k^{2}$ as before and $V=V(r):=\mu^{1/2}rU$ a radial function with $U$ as in (4.13). Then after a straightforward calculation, as a leading order term in a WKB expansion, $V(r)$ must satisfy

[TABLE]

with transmission conditions for $V$ to leading order

[TABLE]

where $\beta^{2}=\rho(r)\mu(r)^{-1}-\omega^{-2}r^{-2}k^{2}$ and $\{r=b\}$ is an interface, 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.

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 $\partial_{r}u=0\text{ on }\partial M$. Indeed, one would merely choose $\mathcal{N}u=\mu\partial_{r}u$ instead without the [math]th order term. However, the boundary condition for $V$ would be of the form

[TABLE]

with $K$ signifying a smooth radial function. Nevertheless, the leading order (in $\omega$) asymptotic behavior for $V$ stays the same despite the $K$ term as clearly seen in the calculation of Appendix B. Thus, our analysis applies with no change using the standard Neumann boundary conditions. This should come as no surprise since in [21], the [math]’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.

In addition, we could also consider a Dirichlet boundary condition, where for $V$, it is also $V=0$ on $\partial M$. This would slightly modify the Debye expansion in Appendix B by constant factors. Nevertheless, the same argument holds to obtain the trace formula and recover the length spectrum. More general boundary conditions such as Robin boundary conditions may be considered as well. However, since we only need to look at the principal term in the high frequency asymptotics, this would just reduce to the Neumann boundary case. Thus, our arguments work with all these boundary conditions, and we choose Neumann boundary conditions only because it has a natural interpretation from geophysics.

An interesting feature of the trace formula in this setup is that a broken ray $\gamma$ can have legs that glide along the interface. This happens when a reflected ray hits an interface at a critical angle leading to a transmitted leg that glides along the interface. Technically, such a ray is not a broken geodesic of the metric $g$, but it will be a limit of periodic broken geodesics as shown in section 2.2 and makes a contribution to the singular support of the trace as an accumulation point.

Since the length spectral rigidity theorems only require the basic length spectrum, the main goal is to determine the leading contribution of basic rays without gliding legs to the trace.

Proposition 4.1**.**

(Non-gliding case) Suppose the radial wave speed $c$ satisfies the extended Herglotz condition and the periodic conjugacy condition (definition 4.7).

Suppose $T=T(p_{\gamma})\in\operatorname{lsp}(c)$ corresponds to a periodic ray $\gamma$ with ray parameter $p_{\gamma}$ such that no periodic ray with a gliding leg has period $T$. Then there exists a neighborhood of $T$ such that, the leading order singularity of $(\operatorname{Tr}(\partial_{t}G))(t)$ near $T(p_{\gamma})$ is the real part of

[TABLE]

where

1. $\bullet$

the sum is taken over all equivalence classes $[\gamma]$ with period $T(p_{\gamma})$ and ray parameter $p_{\gamma}=p_{[\gamma]}$. 2. $\bullet$

$N(p_{\gamma})$* is the Keller-Maslov-Arnold-Hörmander (KMAH) index associated to $\gamma$;* 3. $\bullet$

$c$* independent of $[\gamma]$;* 4. $\bullet$

$\left\lvert SO(3)\right\rvert$* is the volume of the compact Lie group $SO(3)$ under the Haar measure.* 5. $\bullet$

$Q(p_{\gamma})$* is a product of reflection and transmission coefficients of the corresponding broken ray.* 6. $\bullet$

$n(p_{\gamma})\in\mathbb{N}$* is a combinatorial constant counting the number of dynamic analogs of $\gamma$.*

Moreover, if the principal amplitude injectivity condition holds, the distribution $(\operatorname{Tr}\,(\partial_{t}G))(t)=\sum_{n,l}(2l+1)\cos(t{}_{n}\omega_{l})$ is singular at the lengths of periodic basic rays.

Remark 4.2**.**

Our proof will show that one may obtain the leading order contribution of $\gamma^{l}$, which is $\gamma$ traversed $l$ times, from the above expression for $\gamma$. The contribution from $[\gamma^{l}]$ will be

[TABLE]

Remark 4.3**.**

Note the above trace formula is almost identical to that of [14] except for the $Q(p_{\gamma})$ term. This is natural since a wave corresponding to a periodic broken bicharacteristic in this nonsmooth case will have a principal symbol containing transmission and reflection coefficients while the rest of the principal symbol remains the same. The KMAH index also differs slightly than the nonsmooth case when a turning ray grazes an interface.

Remark 4.4**.**

Similar to remark 2.5 in [14], our trace formula holds in an annulus where the boundary is not geodesically convex unlike the case in [21]. Hence, there could be periodic grazing rays at the inner boundary of the annulus or rays that graze an interface. As described in [36], grazing rays are bicharacteristics that intersect the boundary of a layer tangentially, have exactly second order contact with the boundary, and remain in $\bar{M}$. This is another reason our proof is via a careful study of the asymptotics of the eigenfunctions rather than the parametrix construction appearing in [21], where the presence of a periodic grazing ray would make the analysis significantly more technical (cf. [36, 1]). The spherical symmetry essentially allows us to construct a global parametrix (to leading order) to obtain the leading order contribution of a periodic grazing ray to the trace, which would be more challenging in a general setting (see Appendix B and C for the analysis and [6] for a similar computation). The leading order contribution of the grazing ray has the same form as in the above proposition, but the lower order contributions will not have this “classical” form since stationary phase cannot be applied to such terms, and will instead involve Airy functions as in [6] and [14, Appendix B]. Nevertheless, we note that for the main theorems, we do not need to recover the travel time of a periodic grazing ray if one exists. Travel times of sufficiently many non-grazing basic rays suffice. Our methods also produce a precise trace formula where periodic orbits are no longer simple as in [21], but come in higher dimensional families (see [22, 9, 10, 19] for related formulas albeit in different settings).

We showed in section 2.2 that a ray with a gliding leg is a limit of broken non-gliding rays, and we can also describe its contribution to the singular support to leading order. Let $\gamma$ be a periodic broken ray with travel time $T$ and contains a gliding leg (see [8, Figure 4.1] for a diagram of such a ray in the piecewise constant wavespeed setting). By lemma 2.6, there is a sequence of non-degenerate closed broken rays $\gamma_{n}$ with travel time $T_{n}$ such that $T_{n}\nearrow T$ and $\gamma_{n}$ converges to $\gamma$. We will state our trace formula near gliding rays in the same form as [5, Theorem (42)]. Denote $a_{n}=a_{n,[\gamma_{n}]}$ to be the coefficient in (4.20) in front of $(t-T_{n}+i0)^{-5/2}$ corresponding to ray $\gamma_{n}$. We assume that there are no periodic broken rays with travel time $T$ besides for $\gamma$ and its image under the group action. Let us introduce the notation for any real number $s$,

[TABLE]

We will prove the following proposition.

Proposition 4.5**.**

Let $T$ be as above, and let $J$ be a small enough interval containing $T$ such that $\operatorname{lsp}(c)\cap J=\{T_{n}\}_{n=1}^{\infty}\cup\{T\}$.

Then

[TABLE]

where $R(t)$ is a distribution that lies in the Sobolev space $H^{-2+\delta}$ for some $\delta>0$.

Note that this is a genuine error estimate even though we do not have a sharp result on which Sobolev space contains $R(t)$ since the sum in the formula above lies in $H^{-2-}_{loc}$. Proposition 4.5 is not needed for spectral rigidity and will be proved in appendix B.3.1.

Also, implicit in the above proposition is that away from the singularities, the infinite sum converges. It is not clear which Sobolev space $R(t)$ belongs to since we only compute the principal term in the trace (which appears as the sum in the above proposition) using stationary phase, and we show that the remainder is in a weaker Sobolev space even though we cannot use stationary phase for it.

In fact, it is not even clear whether a term of the form $(t-T+i0)^{-\epsilon}$ appears in $R(t)$. Denote $Z(t)=\text{Tr}(\partial_{t}G)(t)$. Then for small enough $\epsilon>0$, $(T-\epsilon,T)\cap\operatorname{lsp}(c)=\{T_{n}\}_{n=1}^{\infty}$ while $(T,T+\epsilon)\cap\operatorname{lsp}(c)=\emptyset$. Thus $\text{Re}Z(t)$ is $C^{\infty}$ for $t\in(T,T+\epsilon)$, and it becomes an interesting question, what is the asymptotic behavior of $Z(t)$ as $t\to T$ from the right? This is subtle and Colin de Verdiére (see [15, 37]) showed how in certain, simpler examples than what we consider here, $Z(t)$ is actually $C^{\infty}$ on $[T,T+\epsilon)$ for some $\epsilon$. Thus, the trace is actually smooth from the right up to and including $T$ (it is obviously not smooth from the left). Červený points out in [8] that the contribution of the singularity precisely at $T$ cannot be investigated with ray theory in this setting, and it remains an open question of the precise nature of this singularity. However, in our computations of the principal term in the WKB expansion, it is not present, which is how we know it can only be in a lower order term, if it is there at all.

The trace formula allows us to recover the basic length spectrum from the spectrum, and then apply the theorems on length spectral rigidity to prove Theorem 1.2.

4.3. Proof of the trace formula

We need several preliminary computations before proving proposition 4.1. The key to the trace formula is the Debye expansion that will give geometric meaning to the leading order amplitudes of the radial eigenfunctions. A key step will be a derivation for an alternative way of expressing $I_{d}$ in (4.7).

4.3.1. A key formula for the Green’s function

As pointed out in [40], the inner product $I_{d}$ can be expressed in terms of the derivatives of a quantity involving the radial eigenfunctions $U_{d}(r)$ as well as their radial derivatives with respect to frequency $\omega$. We repeat the argument here to show that it holds even when the PDE parameters have discontinuities.

The key is obtaining a special formula for $\langle U_{n},U_{n}\rangle$ shown in [35]. We recall the ordinary differential equation (4.3) for the radial constituent of the eigenfunction:

[TABLE]

Here $U=U_{k}=U_{l}$ denotes the above solution for general $\omega$ while $U_{n}$ is a solution for such $\omega_{n}={}_{n}\omega_{l}$ such that $T(U_{n})=\mu(\partial_{r}-r^{-1})U_{n}=0$ at $r=1$ and $r=R$. It will be convenient to write

[TABLE]

Multiply (4.22) by $U_{n}$ and (4.23) by $U$ and subtract the two equation to get

[TABLE]

which may be simplified to

[TABLE]

We integrate over $(R,1)$ to obtain

[TABLE]

Above, we use that $U,U_{n},T,T_{n}$ are continuous across the interface to apply the fundamental theorem of calculus. Let us suppose $\omega$ is not an eigenfrequency and then take the limit as $\omega\to\omega_{n}$. Let

[TABLE]

using the Neumann conditions. Note that the solutions to $D=0$ are precisely the eigenvalues ${}_{n}\omega_{l}$ determined by the Neumann boundary conditions. A key fact is that even for such general solutions, we can enforce the inner boundary condition $T(U)\restriction_{r=R}=0$ to leading order while still keeping $\omega$ generic. This simplifies the computations so that

[TABLE]

Then by L’Hospital’s rule using the limit $\omega\to\omega_{n}$, we obtain

[TABLE]

Next we recall

[TABLE]

Where $I_{d}=I_{n,l}$ is equal to $l(l+1)\int_{r=R}^{1}\rho r^{2}U_{n}^{2}\operatorname{d}\!r.$

What we have shown is that

[TABLE]

so the Green’s function becomes

[TABLE]

Next, observe that ${}_{n}\omega_{l}$ are exactly the zeros of $D$ so we can replace the sum over $n$ by a complex line integral over $\omega$. First use $\text{Re}\tfrac{e^{-i\omega t}}{i}=\sin(\omega t)$. Then for fixed $l$, we compute as in [40]

[TABLE]

where the residue at $\omega={}_{n}\omega_{l}$ of the integrand is calculated via

[TABLE]

and one uses L’Hospital’s rule to get the desired formula. As in [40], the lack of a prefix $n$ on $U_{l}(r)$ and $U_{l}(r^{\prime})$ indicates that these are general solutions which do not necessarily satisfy the free-surface boundary conditions although we are enforcing the inner boundary condition.

Remark 4.6**.**

We note that [14] also used residue theory to compute the infinite sum over $n$. However, the argument would not readily apply here since ${}_{n}\omega_{l}$ is more complicated in our case, so we employ a trick to circumvent using the equations involving ${}_{n}\omega_{l}$, which cannot be solved explicitly.

Thus, we have managed to write $G$ as the Fourier transform in $\omega$ of $D^{-1}\mathbf{D}_{l}(\mathbf{D}_{l})_{0}$. Taking the inverse of the transform, we obtain

[TABLE]

This corresponds with the residue theory in [14] to calculate the infinite series over $n$.

4.3.2. Poisson’s formula for the Green’s function

We abuse notation and denote

[TABLE]

in the formula for $G$ to not treat the curl operations at first. This will not cause risk of confusion since we will specify the exact moment we apply the curl operators. Note that $U_{l}$ does not necessarily satisfy the Neumann boundary conditions.

Proof of proposition 4.1.

By the identical argument in [14, Appendix A], we use Poisson’s formula to rewrite $\hat{G}(x,x_{0},\omega)$ in a different form.

[TABLE]

Note that $H(k)$ has the general eigenfunctions that do not necessarily satisfy Neumann boundary conditions. We substitute $k=\omega p$ so $k\operatorname{d}\!k=p^{-1}\operatorname{d}\!p$ and the above expression becomes (see [14, Appendix A] for details)

[TABLE]

where $Q^{(j)}_{k}(\cos\Theta)$ are certain travelling wave Legendre functions described in [14, Appendix A].

To obtain the leading order asymptotics of the above formula, we will eventually employ the method of steepest descent. Before doing so, we will obtain an expression for $U_{k}(r)$ that has a geometric meaning representing all the multiple scattering of a single ray interacting with not just the boundary, but the interfaces as well. In the Appendix B, we obtain a Debye series expansion of the $D^{-1}H(\omega p)$ term in the above sum.

After a lengthy computation in Appendix B, we write down the updated formula for a single term in the sum over $s$ in the Green’s function from (B.43) when $r$ and $r_{0}$ are in the same layer,

[TABLE]

where the formula is nearly identical to that of [14] with several key differences that encode (using the multiindex $M$) the amplitude and broken ray path consisting of reflecting/transmitting legs. Each component of $M$ indicates the number of reflected or transmitted legs of the ray in a particular layer. First, $Q_{M,i}(p)$ is the leading amplitude of the wave, which is a product of reflection and transmission coefficients corresponding to the legs of a ray connecting two points at $r$ and $r_{0}$ with epicentral distance $\Theta$ (see (B.32) and (B.17)), and $n_{M}$ is a combinatorial coefficient counting the dynamic analogs of this ray. Here, $\tau_{M,i}(r,r_{0};p)$ is the radial component of the travel time of a broken ray with ray parameter $p$ that connect two points at $r$ and $r_{0}$. It is the sum of the radial travel times of each of the reflected and transmitted legs of the ray (see (B.33) and (B.34)). Hence, $\tau_{M,i}$ and $Q_{M,i}$ encode the phase and amplitude (with all the reflections/transmission) of the wave associated to a particular ray. The index $i=1,\dots,4$ corresponds to different ray paths with zero or one reflections connecting the source and receiver located at the radii $r$ and $r_{0}$ analogous to [40]; once we take the trace, and apply the method of steepest descent, only the terms with $i=1,4$ make a contribution to the leading order asymptotics. Moreover, when taking the trace, the terms with $i=1$ and $i=4$ are identical so we will drop the subscript $i$. Also, $N_{M,i}=N_{M,i}(p)$ is the KMAH index which is piecewise constant depending on the value of $p$ and is also affected by a ray grazing an interface.

Method of steepest descent

As in [14, Section 3.2], we carry out the method of steepest descent in the integration over $p$. At this point, the argument is identical so we will be brief. Considering (4.31), we interchange the order of summation and integration, and invoke the method of steepest descent in the variable $p$. Also notice that the path of integration is beneath the real axis, while taking $\omega>0$. We carry out the analysis for a single term, $s=1$. For $s=2,4,\ldots$ we have to add $sp\pi$ to $\tau_{M,i}$, and for $s=3,5,\ldots$ we have to add $(s-1)p\pi$ to $\tau_{M,i}$, in the analysis below.

Considering

[TABLE]

as the phase function (for $s=1$) and $\omega$ as a large parameter, we find (one or more) saddle points for each $i$, where

[TABLE]

Later, we will consider the diagonal, setting $r_{0}=r$ and $\Theta=0$. We label the saddle points by $k$ for each $M,i$ (and $s$). We note that $r,r_{0}$ and $\Theta$ determine the possible values for $p$ (given $M$,$i$ and $s$) which corresponds with the number of rays connecting the receiver point with the source point (allowing conjugate points). Hence, there can be multiple saddle points for a fixed $M,i,s,r,r_{0},\Theta$. For $s=1$, the rays have not completed an orbit. With $s=3$ we begin to include multiple orbits.

We then apply the method of steepest descent to (4.31) with a contour deformation as in [14, Section 3.2] and we obtain

[TABLE]

as $\Theta\to 0$, where

[TABLE]

in which

[TABLE]

The $\tilde{N}_{Mik}$ contribute to the KMAH indices, while the $T_{Mik}$ represent geodesic lengths or travel times. The orientation of the contour (after deformation) in the neighborhood of $p_{k}$ is determined by $\operatorname{sgn}\partial_{p}^{2}\tau_{M,i}\restriction_{p=p_{k}}$. Besides for the geometric spreading factor, the leading order amplitude is $Q_{M,i}(p)$, which is just a product of reflection and transmission coefficients corresponding to the legs of the associated ray; terms involving curvature of the interface do not appear in the lead order term and only make an appearance in the subsequent term that is not necessary for the theorem. We note that

• •

$\tilde{N}_{Mik}=\tilde{N}_{s;Mik}(r,r_{0},\Theta)$ for multi-orbit waves ($s=3,4,\ldots$) includes polar phase shifts produced by any angular passages through $\Theta=0$ or $\Theta=\pi$ as well;

• •

if $r$ lies in a caustic, the asymptotic analysis needs to be adapted in the usual way.

Next, we take the trace of $\partial_{t}G$ by restricting to $(r=r_{0},\Theta=0)$ and integrating. The phase function on the diagonal is $T_{Mik}=\tau_{M,i}(r,r,p_{k})+\pi(s-1)p_{k}$ and we apply stationary phase in the variables $r,\theta,\psi$ with large parameter $\omega$. This is a standard computation exactly as done in [14, Section 3.2].

Following the computation in [14], we obtain the leading order term in the trace formula as

[TABLE]

where $L_{k}$ is the travel time of a ray with only transmitted legs from $r=1$ to $r=R^{*}$ (see (2.2)). Note that the critical set becomes $\Theta_{M,k}=\partial_{p}\tau_{M}(p_{k})$ so $\partial^{2}_{p}\tau_{M}(p_{k})=\partial_{p}\alpha_{M,k}$ when restricting to the diagonal. Also, we use that here,

[TABLE]

is independent of $r$ on the critical set. We note that $p_{k}$ exists only for $\left\lvert M\right\rvert$, and $s$, sufficiently large, which reflects the geometrical quantization.

∎

Harmonics of the principal ray

From the argument above, if $\gamma$ is a periodic orbit with period $T_{s,Mik}$ for some indices $s,M,i,k$ described above, the principal symbol of the contribution of $\gamma$ to the trace is as above. We can immediately write down the leading order contribution of $\gamma^{l}$ which is $\gamma$ travelled $l$ times. The travel time will be $lT_{s,Mik}$ . Then $Q_{M,i}(p_{k})$ becomes $Q_{M,i}(p_{k})^{l}$, $M_{ik}$ becomes $lM_{ik}$, and $p_{k}^{-2}\partial_{p}\alpha_{M,ik}$ becomes $lp_{k}^{-2}\partial_{p}\alpha_{M,ik}$.

4.4. Spheroidal modes

The above trace formula, theorems 1.6 and 1.2 are essentially dealing with a scalar wave equation with a single wavespeed. The analysis for toroidal modes reduced to a scalar wave equation with an associated Laplace-Beltrami operator. However, our methods can also treat the isotropic elastic setting to include spheroidal modes (with the PDE described in [13, Chapter 8]) where two wavespeeds ($c_{P}$ and $c_{S}$) are present corresponding to the $P$-waves and the $S$-waves, and there is a spectrum associated to the elliptic, isotropic elastic operator. In the elastic setting, each leg of a broken geodesic will be a geodesic for either the metric $c_{P}^{-2}dx^{2}$ or $c_{S}^{-2}dx^{2}$, so there is an associated length spectrum as well that includes mode converted legs. Thus, theorem 1.6 can be extended to the case of the elastic operator by using corollary 1.4 if the length spectrum (or a dense subset) can be recovered by a trace formula from the spectrum. The theorem would take the form

Theorem 4.7** (Elastic spectral rigidity with moving interfaces).**

Fix any $\varepsilon>0$ and $K\in\mathbb{N}$, and let $c_{P,\tau}(r)$ and $c_{S,\tau}(r)$ be an admissible family of profiles with discontinuities at $r_{k}(\tau)$ for all $k=1,\dots,K$. Suppose that the length spectrum for each $c_{P/S,\tau}$ is countable in the ball $\bar{B}(0,1)\subset\mathbb{R}^{3}$. Assume also that the length spectrum satisfies the principal amplitude injectivity condition and the periodic conjugacy condition.

Suppose $\operatorname{spec}(\tau)=\operatorname{spec}(0)$ for all $\tau\in(-\varepsilon,\varepsilon)$. Then $c_{P,\tau}=c_{P,0}$, $c_{S,\tau}=c_{S,0}$ and $r_{k}(\tau)=r_{k}(0)$ for all $\tau\in(-\varepsilon,\varepsilon)$ and $k=1,\dots,K$.

Thus, all we need is to extend proposition 4.1 to the elastic case and then apply corollary 1.4. Since the calculation is similar, but a more tedious version of the case we consider here, We will just provide an outline of the proof.

1. (1)

The Green’s function associated to just the spheroidal modes can be computed analogously as in [40, Equation (31)]. 2. (2)

One can then obtain (vector-valued) WKB solutions to approximate spheroidal modes, which are eigenfunctions of the static, elastic operator as in [40, Appendix A] and [13, Chapter 8]. 3. (3)

We can use the methods presented here (with the method of steepest descent for the asymptotic analysis) to then determine the leading order asymptotics of the sum of eigenfunctions to obtain a corresponding trace formula. The scattering coefficients will be determined by the elastic transmission condition, with an associated Debye expansion as done in appendix B. Afterward, the stationary phase analysis will lead to the same form as (4.20) but the reflection and transmission coefficients appearing in $Q(p_{\gamma})$ will be different to account for mode conversions. Also, $\alpha(p_{\gamma})$ will be modified with the appropriate wave speed appearing in each constituent of the linear combination of epicentral distances that correspond to an associated $P$ or $S$ leg of $\gamma$. 4. (4)

The computation in [40] does not treat glancing nor grazing rays, but their formulas can be modified with the methods presented here to account for such rays as well. The $n(p_{\gamma})$ appearing in (4.20) will again count the number of “dynamic analogs” associated to $\gamma$ as described in [24] for the spheroidal case; that paper also has several figures of broken geodesics in the spheroidal case.

Under an analog to the principal amplitude injectivity condition for spheroidal modes, one can recover the basic length spectrum for each of the two wave speeds. One then uses Corollary 1.4 to recover both wave speeds.

5. Declarations

Funding

MVdH was supported by the Simons Foundation under the MATH + X program, the National Science Foundation under grant DMS-1815143, and the corporate members of the Geo-Mathematical Imaging Group at Rice University. JI was supported by the Academy of Finland (projects 332890 and 336254).

Conflict of interest/Competing interests

Financial interests: The authors declare they have no financial interests.

Non-financial interests: The authors declare they have no non-financial interests.

Availability of data and material

Not applicable

Code availability

Not applicable

Appendix A Class of metrics satisfying (A1)-(A4)

Here, we present a class of metrics that meet all the assumptions mentioned in the theorems. Demonstrating the satisfaction of assumption (A.2) on a layered manifold (or any analogous concept of a simple length spectrum) is a challenging task, even for the most basic Herglotz metrics. Therefore, we instead prove that this assumption is generic within our class of metrics. We subsequently prove that the Euclidean ball possesses a simple length spectrum, thereby establishing the non-emptiness of the class of metrics satisfying our assumptions.

It is worth noting that the proof of Theorem 1.6 only necessitates the recovery of the basic length spectrum, and one may omit any finite number of such lengths. Thus, we can allow for cancellations in the trace formula as long as they do not occur infinitely many times for periods in the basic length spectrum.

We will begin by examining a single interface at $r=\mathbf{f}$ and then proceed to generalize this approach to additional layers. Let us select a function $c(r)$ such that

[TABLE]

and

[TABLE]

where $a,a^{\prime},b,b^{\prime}>0$ are parameters, ensuring that condition (A4) is satisfied. Since the proof of the main theorem relies solely on the basic length spectrum, we will verify our assumptions for this specific case.

For basic rays in the upper layer, we can calculate the epicentral distance of a single leg as follows:

[TABLE]

Here, the value of $R^{*}$ depends on the parameter $p$. We can further break down the calculation:

[TABLE]

where in the case of reflecting rays, $R^{*}=\mathbf{f}$ in the first layer, and for diving rays (as discussed in Appendix B), $R^{*}$ is determined by the condition $(c(R^{*})/R^{*})^{-2}-p^{2}=0$. In this scenario, we have:

[TABLE]

This expression only equals zero when $p^{2}=a/2$. We can similarly verify that $\alpha^{\prime}(p)$ does not vanish in the reflecting region ($0<p<\sqrt{a+b\ln(\mathbf{f})}$) except for $p^{2}=a(a+b\ln(\mathbf{f}))/(2a+b\ln(\mathbf{f}))$. A similar argument holds for basic rays in the second layer, and this argument readily extends to any finite number of layers, each having a metric within the above class. If we have $N$ layers with similar such metrics, there is an associated $\alpha_{i}(p)$ for layer $i$. Then for fixed nonnegative integers $n_{1},\dots n_{N}$, the epicentral distance is $\alpha(p)=\sum n_{j}\alpha_{j}$ so $\alpha^{\prime}(p)=0$ for at most finitely many values of $p$. Thus, (A3) is satisfied

To get a periodic basic ray, we need $\alpha=2\pi q$ for some $q\in\mathbb{Q}$. For diving rays, this implies

[TABLE]

Solving for $p$, we obtain

[TABLE]

Hence, periodic diving rays are in one-to-one correspondence with elements of $\mathbb{Q}$ such that

[TABLE]

If we choose $a,b$ such that $a/b\notin\pi\mathbb{Q}$, then $p^{2}\neq a/2$ for a periodic ray, and the clean intersection hypothesis is satisfied on all diving periods. A similar calculation provides the result for reflecting rays and rays in the second layer. Periodic reflecting rays in a layer are in one-to-one correspondence with elements of $\mathbb{Q}$ such that $ad-\pi^{2}q^{2}b^{2}\geq 0$ where $d=a+b\ln\mathbf{f}$. If we restrict to $a,b$ such that $\sqrt{ad}/b\notin\pi\mathbb{Q}$, then the clean intersection hypothesis is satisfied for all basic reflecting rays. The extension to any finite number of layers follows analogously since we only require the basic length spectrum.

It is important to note that this argument can also extend to non-basic rays since assumption (A.1) could fail on at most finitely many primitive periodic rays, which still suffices for length spectral rigidity.

Remark A.1**.**

It should be noted that the proof of the main theorems doesn’t require the recovery of all lengths in $\operatorname{lsp}(c)$. Therefore, even if the clean intersection hypothesis fails for a finite number of rays and those lengths are not recovered with the trace formula, the theorems remain valid with the same proof.

A.1. Genericity of assumption (A.2)

To meet the requirements of (A.2), we initially consider two layers with a wave speed given having the form $c=a+b\ln(r)$ in each layer, and prove that assumption (A.2) is generic (in a sense described below) among metrics in this class. The argument we provide will readily extend to the case of a finite number of layers. Let $M$ represent the annulus where we connect both of these layers (layer $1$ and layer $2$), with wave speed $c_{1}$ in the upper layer and $c_{2}$ in the lower layer. Let $c$ denote the piecewise smooth wave speed with the length spectrum $\operatorname{lsp}$.

Suppose we have a closed basic ray $\gamma$ with a period $T=T(p_{\gamma})$ such that a cancellation occurs in the trace formula. Let’s outline the conditions under which this is possible. First, denote the closed ray $\beta=\beta(p)$ that corresponds to such a cancellation in the trace formula where $[\gamma]\neq[\beta]$. Hence, the principal term in the wave trace corresponding to $[\gamma]$ cancels with that of $[\beta]$. Let $L_{i}=L_{i}(p)$ denote the length of one leg in layer $i$ for ray $\beta$. There exist $n_{1}$ and $n_{2}$ nonnegative integers that satisfies the equation:

[TABLE]

Furthermore, to ensure periodicity, $n_{1},n_{2},p$ must also satisfy:

[TABLE]

for some integer $m$.

Definition A.2**.**

We will define a period $T=T(p_{\gamma})$ as not simple if condition (A.1) is satisfied for some $n_{1},n_{2},q,m,$ and closed ray $\beta$ with a ray parameter $p$ such that $[\beta]\neq[\gamma]$, as described above. These conditions represent the simultaneous requirements necessary for a cancellation in the wave trace. In other words, if a period $T$ is simple, then it will not result in a cancellation in the wave trace. We will say that $\operatorname{blsp}(c)$ is simple if each non-gliding period in this set is simple. A period is gliding when its associated ray has a gliding leg, and such periods are not needed to prove length spectral rigidity.

Suppose we have a class of metrics

[TABLE]

Our goal is to demonstrate that

[TABLE]

is generic in the Baire category sense. In other words, within the class $\mathcal{C}$, metrics generically satisfy (A1)-(A4). Accumulation points in the length spectrum can pose difficulties, so it’s helpful to define a set for $\delta>0$ as follows:

[TABLE]

Note that the gliding ray traveling once around the outer boundary or interface has a length of $2\pi$ or $2\pi\mathbf{f}$, respectively. Therefore, $\mathcal{G}_{\delta}$ represents lengths that are at least $\delta$ away from each gliding ray along the outer boundary or interface. Thus, throughout the arguments below, closed rays will have periods that stay a fixed distance away from any period associated to a gliding ray.

Let us define

[TABLE]

so that $E=\cap_{n}E_{n}$. Then we have

[TABLE]

We would like to prove that $E_{n}$ is open and that $E_{n}$ is dense in $E_{n-1}$ for each $n$ to apply the Baire category theorem. This will suffice since the proof of the main theorems go through even if we do not recover periodic lengths in the lengths spectrum that align with a gliding period. Furthermore, we demonstrated above that within this class of metrics, only a finite number of degenerate periods can exist. The spectral rigidity proof holds even if we do not recover any finite subset of the basic length spectrum. Therefore, all metrics in the set $E$ satisfy the assumptions necessary for theorem 1.6 to hold.

We will follow similar arguments to [2]. Using the notation there, for a metric $g$, denote $\{\Phi_{t}^{g}\}$ to be the bicharacteristic flow on $T^{*}M$, which restricts to a bicharacteristic flow on $S^{*}M.$ For a particular ray $\gamma$, we denote $N_{refl}=N_{refl}([\gamma])=(N^{(1)}_{refl},N^{(2)}_{refl}),$ where

[TABLE]

and

[TABLE]

A.1.1. Openness

To prove that $E_{k}$ is open, let $g_{c}=g_{c_{\mathbf{a}}}\in E_{k}$ with $\mathbf{a}=(a,b,a^{\prime},b^{\prime})$. We want to demonstrate that there is a neighborhood $O\subset\mathbb{R}^{4}$ of $\mathbf{a}$ such that $g_{c_{\mathbf{b}}}\in E_{k}$ for each $\mathbf{b}\in O$.

Assume, for the sake of contradiction, that this is not the case. Then there exists a sequence $\mathbf{a}_{n}\to\mathbf{a}$ and a sequence of periods $T_{n}=T_{n}(p_{n})\in\operatorname{blsp}(c_{\mathbf{a}_{n}})$, each of which is not simple and has a length $\leq k$ and belongs to $\mathcal{G}_{1/k}$. We will show that a subsequence of these periods converges to a period in $\operatorname{blsp}(c_{\mathbf{a}})$ that is also not simple, leading to a contradiction.

By definition, for each $n$, there exists closed rays $\gamma_{n},\tilde{\gamma}_{n}$ with associated $p_{n},\tilde{p}_{n}$ such that $[\gamma_{n}]\neq[\tilde{\gamma}_{n}]$, but the associated period $\tilde{T}_{n}=T(\tilde{p}_{n})\in\operatorname{lsp}(c_{n})$ is such that:

[TABLE]

Each $p_{n}$ is associated to a periodic ray equivalence class, so there is a $v_{n}\in S^{*}M$ such that

[TABLE]

Likewise, we have $\tilde{v}_{n}\in S^{*}M$ such that

[TABLE]

By compactness and passing to subsequences (while still denoting the index by $n$), we can find $T,\tilde{T}\leq k$, $v,\tilde{v},p,\tilde{p}$ such that $T_{n}\to T$, $\tilde{T}_{n}\to\tilde{T}$, $p_{n}\to p$, $\tilde{p}_{n}\to\tilde{p}$, $v_{n}\to v$, and $\tilde{v}_{n}\to\tilde{v}$. Since we know $g_{n}\to g_{c}$, upon taking limits in the equations above, we obtain:

[TABLE]

and

[TABLE]

and $T(p)=\tilde{T}(\tilde{p})$. Let $\gamma,\tilde{\gamma}$ denote the rays $\Phi^{g_{c}}_{t}v$ and $\Phi^{g_{c}}_{t}\tilde{v}$ projected to $M$ respectively. They are periodic with periods $T$ and $\tilde{T}$, respectively. The assumption regarding $\mathcal{G}_{1/k}$ implies there exists $\delta>0$ such that $L_{1}(\tilde{p}_{n})>\delta$ and $L_{2}(\tilde{p}_{n})>\delta$ for large enough $n$ and assuming the ray is in both layers. If the ray is basic, then a similar argument will hold. By (A.1), we have

[TABLE]

for some nonnegative integers $z_{1,n},z_{2,n}$. This shows that $N_{n,refl},N_{n,WN}$ remain bounded for large enough $n$, and in particular, become constant. Thus, $N_{n,refl}=N_{refl}([\gamma])$ and $N_{n,WN}=N_{WN}([\gamma])$ for large enough $n$. Let us show that $[\gamma]\neq[\tilde{\gamma}]$. For the sake of contradiction, suppose $[\gamma]=[\tilde{\gamma}]$. This would imply that $N_{refl}=\tilde{N}_{refl}$ and $N_{WN}=\tilde{N}_{WN}$. Hence, $N_{n,refl}=\tilde{N}_{n,refl}$ and $N_{n,WN}=\tilde{N}_{n,WN}$ for large enough $n$. Since $T_{n}=\tilde{T}_{n}$, then $\gamma$ and $\tilde{\gamma}$ are kinematic analogs so that $[\gamma_{n}]=[\tilde{\gamma}_{n}]$. This is because since $p=\tilde{p}$, then $p_{n}$ is near $\tilde{p}_{n}$ for large $n$, and due to nondegeneracy, if $T(\sigma)$ denotes the travel time for $\sigma$ near some $p_{n}$, then $\partial_{\sigma}\restriction_{\sigma=p_{n}}T\neq 0$. This implies $T(p_{n})\neq T(\tilde{p}_{n})$ when $p_{n}\neq\tilde{p}_{n}$ for large enough $n$. Thus, for large $n$, $[\gamma_{n}]=[\tilde{\gamma}_{n}]$, which is a contradiction. Hence, $[\gamma]\neq[\tilde{\gamma}]$ so $T$ is a non-simple period for $g_{c}$, which is a contradiction. This establishes openness.

A.1.2. $E_{k+1}$ is dense in $E_{k}$

Let $g_{c}\in E_{k}$. We will show there exists a $\tilde{g}\in E_{k+1}$ arbitrarily close to $g_{c}$.

First, let us assume that $M$ is a single layer, and describe afterward how the proof extends to multiple layers. Let $[\gamma_{1}],\dots[\gamma_{N}]$ be all the distinct closed rays equivalence classes whose periods are $\leq k+1$, belong to $\mathcal{G}_{k+1}$, and are non-degenerate. In view of the non-degeneracy computation at the start of the section, and that we are a fixed distance $1/(k+1)$ away from the periods of gliding rays, the proof in [2] shows that there are finitely many such trajectories.

For each $i$, we have an associated epicentral distance $\alpha(p_{i})$ for the ray $\gamma_{i}$. Consider any analytic perturbation of the metric $g_{c}$ denoted $g_{c_{h}}$ with parameter $h$. By [14], we may find a $\epsilon>0$ and a function

[TABLE]

such that the ray

[TABLE]

has associated

[TABLE]

and is thus periodic. Also $N_{i,refl,h}$ and $N_{i,WN,h}$ remain constant as well under varying $h$.

By the proof in [2], for any period ray $\gamma_{h}$ associated to metric $g_{c_{h}}$ whose period is $\leq k+1$, is non-degenerate, and belongs to $\mathcal{G}_{k+1}$, then $[\gamma_{h}]$ must coincide with one of

[TABLE]

Hence, if $g_{c}$ is not in $E_{k+1}$, then there exists $T,\tilde{T}\in\operatorname{blsp}(c)$, such that $T(p)=\tilde{T}(\tilde{p})$ with length $\leq k+1$, they are non-degenerate, and belong to $\mathcal{G}_{k+1}$. We will show that under a slight perturbation of $c$, $T(p)$ becomes simple while keeping all the simple periods of $g_{c}$ simple under perturbation. Without loss of generality, we may assume these two periods correspond to the ray $\gamma_{1}$ and $\gamma_{2}$, and the other periods $T_{3},T_{4}\dots,T_{N}$ are simple (if other periods were non-simple, we can just repeat the following argument until all are simple).

Under a perturbation of $c$, denoted $c_{h}$ with parameter $h$, the perturbed periods $T_{3,h},\dots,T_{N,h}$ remain simple for $h$ small enough due to continuity of $h\mapsto T_{i,h}$ (in fact, differentiability by [14]), the openness of $E_{k}$, and since these are all the periods $\leq k+1$.

Thus, we need to show that we may pick a particular perturbation so that $T_{1,h}$ becomes simple and the proof of density is complete. The perturbed travel time is

[TABLE]

and $\tilde{T}_{h}=T_{2}(p)+h\delta T(\tilde{p})+O(h^{2})$.

If we can arrange that $\delta T(p)\neq\delta T(\tilde{p})$, then $T_{h}\neq\tilde{T}_{h}$ for small enough $h$. Note that for a fixed $\mathbf{a}^{\prime}\in\mathbb{R}^{4}$, consider the analytic perturbation $\mathbf{a}_{h}=\mathbf{a}+h\mathbf{a}^{\prime}$. By openness, $c_{h}\in E_{k}$ for small enough $h$. Since we have an explicit formula for $T_{h}$ and $\tilde{T}_{h}$, we can compute explicitly $\frac{d}{dh}\restriction_{h=0}T_{h}$ and so too for $\tilde{T}_{h}$, and we can pick $\mathbf{a}^{\prime}$ to guarantee that $\delta T$ and $\delta\tilde{T}$ are different. Thus $T_{h}(p_{h})\neq\tilde{T}_{h}(\tilde{p}_{h})$. This establishes the required density.

The above argument readily applies to the case of multiple interfaces since it only relies on having a discrete length spectrum away from a fixed distance away from the gliding periods, and that the quantities $N_{n,\text{refl}}$ and $N_{n,\text{WN}}$ become constant for large $n$ when a family of periodic rays converges to another ray. Thus, the only issue is nonbasic periodic rays that have a gliding leg. However, if we restrict our class of metrics such that

[TABLE]

then such nonbasic rays with a gliding leg do not exist. Hence, the only accumulation points in the length spectrum lie in $\mathcal{G}_{\delta}$, and we do not need to recover such periods for the proof.

On the other hand, if one wanted to consider such nonbasic rays with one or more gliding legs, then these would not create a cancellation in the wave trace, even if their period coincided with the period of a ray in a disjoint equivalence class. This is because the proof of the trace formula shows that the term in the wave trace corresponding to such a ray is smoother than the terms corresponding to non-gliding rays.

A.2. Simplicity of length spectrum of a Euclidean ball

We give a concrete example of a simple basic length spectrum, the Euclidean disc. The basic length spectrum is simple and the full length spectrum is almost simple; the only coincidence is that going twice around a hexagon takes as long as going three times back and forth along a diameter.

The proof relies heavily on number theory. The nature of the proof in the geometrically simplest case suggests that finding any concrete examples in non-Euclidean geometry will be very difficult. This is why genericity arguments like the one presented above are necessary.

Theorem A.3**.**

Let $\gamma_{1}$ and $\gamma_{2}$ be two periodic broken rays in the Euclidean unit disc. Suppose that the rays are basic and not rotations of each other. Then their lengths are different.

In fact, the ratio of their lengths is rational if and only if $\gamma_{1}$ is a hexagon (length $6$) and $\gamma_{2}$ is a diameter (length $4$, counted back and forth).

The proof is easy once we establish the following proposition.

Proposition A.4**.**

If $r_{i}$ for $i=1,2$ are two different rational numbers in $(0,\frac{1}{4}]$, then

[TABLE]

unless $r_{1}=\frac{1}{4}$ and $r_{2}=\frac{1}{12}$ or vice versa.

Proof of theorem A.3.

Each $\gamma_{i}$ is composed of line segments with an opening angle $\alpha_{i}$ as seen from the origin. The condition for periodicity is that $\alpha_{i}\in\pi\mathbb{Q}$, and in Euclidean geometry $\alpha\in(0,\pi]$. We thus have $\alpha_{i}=2\pi\cdot\hat{r}_{i}$ with $\hat{r}_{i}=p_{i}/q_{i}\in(0,\frac{1}{2}]\cap\mathbb{Q}$. The length of one line segment is $l_{i}=2\sin(\alpha_{i}/2)=2\sin(\hat{r}_{i}\pi)$ and the length of the whole ray is $L_{i}=q_{i}l_{i}$.

The claim is that $q_{1}l_{1}\neq q_{2}l_{2}$. In most cases this follows from proposition A.4. The only case not covered by it is when $\hat{r}_{1}=\frac{1}{2}$ and $\hat{r}_{2}=\frac{1}{6}$ or vice versa. These correspond to a diameter (with $q_{1}=2$ and $L_{1}=4$) and a hexagon (with $q_{1}=6$ and $L_{1}=6$). In this one case the ratio is rational, but still $L_{1}\neq L_{2}$, completing the proof of the main claim and the additional claim. ∎

To prove the proposition we will use some notation and concepts from number theory without explanation in this section to keep exposition concise. In what follows, $\zeta_{k}=e^{2\pi i/k}$ is a $k$th root of unity, $\mathbb{Q}(\zeta_{k})$ is the $k$th cyclotomic field, $\varphi$ is Euler’s totient function, $\gcd(a,b)$ is the greatest common divisor of $a$ and $b$, and the degree of a field extension of an element thereof is the degree of its minimal polynomial over the subfield.

Lemma A.5**.**

Let $a$ and $b$ be integers with $a>1$ and $b>1$. Then $\varphi(ab)=\varphi(b)$ if and only if $b$ is odd and $a=2$.

Proof.

Let $p_{1},\dots,p_{m}$ be the primes appearing in the prime factorizations of $a$ and $b$, so that

[TABLE]

where some of the powers $k_{i}$ and $l_{i}$ may be zero.

If we define the functions $f_{i}\colon\mathbb{N}\to\mathbb{N}$ by

[TABLE]

then

[TABLE]

Since each $f_{i}$ is increasing, we have $f_{i}(l_{i})\leq f_{i}(k_{i}+l_{i})$ for all $i$ and so $\varphi(b)\leq\varphi(ab)$.

The equality $\varphi(b)=\varphi(ab)$ is only possible when $f_{i}(l_{i})=f_{i}(k_{i}+l_{i})$ for all $i$. A quick inspection reveals that $f_{i}$ is strictly increasing with one exception: $f_{i}(0)=f_{i}(1)$ if $p_{i}=2$.

We assumed $a>1$, so at least one $k_{i}$ is non-zero. But $f_{i}(l_{i})=f_{i}(k_{i}+l_{i})$ for $k_{i}>0$ is only possible when $l_{i}=0$ and $k_{i}=1$ and $p_{i}=2$. Therefore $\varphi(b)=\varphi(ab)$ implies that $a=2$ and $b$ is odd. ∎

Lemma A.6**.**

Let $a$ and $b$ be positive integers. The following are equivalent:

1. (1)

$\mathbb{R}\cap\mathbb{Q}(\zeta_{a})=\mathbb{R}\cap\mathbb{Q}(\zeta_{b})$. 2. (2)

Any of the following holds:

1. (a)

$a=b$. 2. (b)

$a=2b$* and $b$ is odd.* 3. (c)

$b=2a$* and $a$ is odd.* 4. (d)

$a$* and $b$ are both in the set $\{1,2,3,4,6\}$.*

Proof.

We only prove “$\implies$” — the reverse implication “$\impliedby$” is simpler and can be easily constructed from the arguments presented here. Suppose that $\mathbb{R}\cap\mathbb{Q}(\zeta_{a})=\mathbb{R}\cap\mathbb{Q}(\zeta_{b})$ and, without loss of generality, that $a<b$.

Intersections of cyclotomic fields satisfy

[TABLE]

so the assumption implies also that $\mathbb{R}\cap\mathbb{Q}(\zeta_{a})=\mathbb{R}\cap\mathbb{Q}(\zeta_{b})=\mathbb{R}\cap\mathbb{Q}(\zeta_{g})$, where we have denoted $g=\gcd(a,b)$ for short. The assumption $a<b$ implies $g<b$.

Since $\zeta_{1}=1$ and $\zeta_{2}=-1$, we have simply $\mathbb{Q}(\zeta_{1})=\mathbb{Q}(\zeta_{2})=\mathbb{Q}$. When $n\geq 3$, the number $\zeta_{n}$ is not real and its degree is $\varphi(n)$ and the minimal polynomial is the $n$th cyclotomic polynomial. The degree of the real number $z_{n}\coloneqq\frac{1}{2}(\zeta_{n}+\overline{\zeta_{n}})$ is $\varphi(n)/2$ as shown in [38]. Because $\mathbb{R}\cap\mathbb{Q}(\zeta_{n})=\mathbb{Q}(z_{n})$, our assumption implies that the degrees of $z_{a}$ and $z_{b}$ and $z_{g}$ agree.

Let us split the analysis in different cases, using the fact that $\varphi^{-1}(1)=\{1,2\}$ and $\varphi^{-1}(2)=\{3,4,6\}$:

1. (1)

$b=5$ or $b\geq 7$: In this case $\varphi(b)>2$ and so $\deg(z_{b})>1$. Thus also $\deg(z_{g})>1$, which implies that $g=5$ or $g\geq 7$. Therefore the formula $\deg(z_{n})=\varphi(n)/2$ is valid for $n=b,g$ and we find that $\varphi(g)=\varphi(b)$. Now that $g|b$ and $g<b$, lemma A.5 gives that $b=2g$ and $b$ is odd. The only possible value for $a$ is $g$. 2. (2)

$b\in\{2,3,4,6\}$: In this case $\deg(z_{b})=1$ and so $\mathbb{R}\cap\mathbb{Q}(\zeta_{b})=\mathbb{Q}$. The set of values of $a$ for which $\deg(z_{a})=1$ is $\{1,2,3,4,6\}$.

This is exactly the claim. ∎

Proof of Proposition A.4.

We switch from the parameter $r\in(0,\frac{1}{4}]$ to $s=\frac{1}{4}-r\in[0,\frac{1}{4})$ so that $\sin(2r\pi)=\cos(2s\pi)\eqqcolon c_{s}$.

Suppose $c_{s_{1}}/c_{s_{2}}\in\mathbb{Q}$. This implies that $c_{s_{1}}\in\mathbb{Q}(c_{s_{2}})$ and so $\mathbb{Q}(c_{s_{1}})\subset\mathbb{Q}(c_{s_{2}})$. The same argument works in reverse, so $\mathbb{Q}(c_{s_{1}})=\mathbb{Q}(c_{s_{2}})$.

Consider first the case when $s_{i}\neq 0$. If $s_{i}=p_{i}/q_{i}$ with $\gcd(p_{1},q_{i})=1$, let us denote $\hat{s}_{i}\coloneqq 1/q_{i}$. The assumption $s_{i}<\frac{1}{4}$ implies that $q_{i}\geq 5$.

By [38] the algebraic numbers $c_{s_{i}}$ and $c_{\hat{s}_{i}}$ have the same minimal polynomial. Therefore

[TABLE]

By lemma A.6 this implies that one of the following holds:

1. (1)

$q_{1}=q_{2}$. 2. (2)

$q_{1}=2q_{2}$ and $q_{2}$ is odd. 3. (3)

$q_{2}=2q_{1}$ and $q_{1}$ is odd.

The other cases of the lemma are ruled out by the property $q_{i}\geq 5$, and the third case listed above can be ignored by symmetry.

Let $P_{q}$ be the set of numbers $p$ with $1\leq p<q/2$ and $\gcd(p,q)=1$. We have $\#P_{q}=\varphi(q)/2$. The cosines $c_{p/q}$ with $p\in P_{q}$ share the same minimal polynomial and they form a basis for $\mathbb{Q}(c_{1/q})$ over $\mathbb{Q}$. They are therefore linearly independent over $\mathbb{Q}$.

Consider first the case $q_{1}=q_{2}\eqqcolon q$. The numerators $p_{1}$ and $p_{2}$ are distinct elements of $P_{q}$ by assumption. The linear independence of the cosines contradicts $c_{s_{1}}/c_{s_{2}}\in\mathbb{Q}$.

Consider then the case $q_{1}=2q_{2}\eqqcolon 2q$ with odd $q$. Because $p_{1}$ and $q$ are odd, there is an integer $\tilde{p}_{1}$ so that $p_{1}=q_{2}-2\tilde{p}_{1}$. It satisfies $\tilde{p}_{1}>q/4>p_{2}$. The trigonometric identity $\cos(x)=-\cos(\pi-x)$ gives $c_{p_{1}/q_{1}}=c_{(q_{2}-2\tilde{p}_{1})/2q_{2}}=-c_{\tilde{p}_{1}/q_{2}}$. Because $\tilde{p}_{1}\neq p_{2}$, the argument of the previous case shows that the ratio of the cosines is irrational, a contradiction.

The only remaining case is $s_{i}=0$. Without loss of generality, take $s_{1}=p/q>0$ (with $\gcd(p,q)=1$ and $q\geq 5$) and $s_{2}=0$. We have $c_{0}=1$, so we are led to study when $c_{p/q}\in\mathbb{Q}$. Applying lemma A.6 in the case $a=1$ and $b=q$ gives that $\mathbb{Q}(c_{p/q})=\mathbb{Q}$ if and only if $q\in\{1,2,3,4,6\}$. Due to this and the bound $q\geq 5$ the assumption $c_{s_{1}}/c_{s_{2}}\in\mathbb{Q}$ gives $q=6$. The only possible value left for $p$ is $1$.

We have found that $s_{1}=\frac{1}{6}$ and $s_{2}=0$ or vice versa. This means that $r_{1}=\frac{1}{4}-\frac{1}{6}=\frac{1}{12}$ and $r_{2}=\frac{1}{4}-0=\frac{1}{4}$ or vice versa, as claimed. ∎

Appendix B Generalized Debye Expansion

In this appendix, we will reduce equation (4.30) into the form (B.43), which resembles a wave propagator. We evaluate $D^{-1}U_{l}(r)U_{l}(r_{0})$ appearing in (4.30) in such a way as to relate it to a certain wave propagator analogous to the computation in [14]. However, the methodology will be different and more laborious to account for the multiple scattering created by the interfaces.

B.1. Single interface case

For simplicity, we first consider a $2$-layered sphere with an upper layer $\Omega_{+}$ and a lower layer $\Omega_{-}$. The general case will follow easily by a recursive argument. The wave speed and density in a layer $\pm$ (region $\pm$) is $c_{\pm},\rho_{\pm}$. We have the upper surface layer $r=1$ the inner boundary $r=R$ (where $R=0$ if we consider the case of a ball) and the interface $r=b$.

Suppose $h^{(1)}(r),h^{(2)}(r)$ are two linearly independent solutions to the second order ODE (4.3) not necessarily satisfying any boundary condition. They implicitly depend on $k$ and $\omega$. Suppose $r>b$ is the $+$ region and $r<b$ is the $-$ region. Write the solutions in the $\pm$ region

[TABLE]

For $u_{-}$, we think of $C$ as being determined by either an inner boundary condition or another transmission condition if there were more layers, and not by the transmission conditions at $r=b$, which instead determine $B$. Hence, to emphasize this point, and to make the computation cleaner, denote

[TABLE]

It will also be useful to consider the solutions to the simpler ODE for $V=\mu^{1/2}rU$ with Neumann boundary condition $\partial_{r}V\restriction_{r=1,R}=0$.

We will think of $h_{+}^{(2)}$ as an incoming wave into the interface $r=b$ and $h_{+}^{(1)}$ as a scattered wave even though the notation is merely symbolic at this point. Without writing it explicitly, the constants are not functions of $r$ but they do depend on $l$ and $\omega$. When we make the substitution, $k=\omega p$, then we will have $A=A(\omega,p)$, and for the remaining constants as well. The general interface conditions to leading order asymptotics as $\omega\to\infty$ are (with $d_{\pm}$ some parameter on each side of the interface; in our particular case, $d_{\pm}=\mu_{\pm}(b)$)

[TABLE]

To ease notation, omit the evaluation at $r=b$ so we have

[TABLE]

Then

[TABLE]

Thus,

[TABLE]

Factor out $jd_{-}h^{(2)}_{+}$ in the numerator and $jd_{-}h^{(1)}_{+}$ in the denominator to get

[TABLE]

Let us use the notation

[TABLE]

We then have

[TABLE]

Following [29, Appendix], we solve for reflection and transmission coefficients in terms of the above functions. We write

[TABLE]

and

[TABLE]

Here, $u_{+}$ and $u_{-}$ are solutions unrelated to the previous $u_{\pm}$. One should think of them as being defined only locally near an interface. The notation is that $R_{++}$ is reflection from above the interface and $T_{+-}$ is transmission from the upper layer $(+)$ to the lower layer $(-)$. The transmission conditions give

[TABLE]

Thus

[TABLE]

So

[TABLE]

So

[TABLE]

Likewise, we can show that

[TABLE]

Debye expansion for A in (B.4)

Using (B.4) and the formulas for reflection and transmission coefficients, we follow Nussenszweig [29](all functions are evaluated at $r=b$ without explicitly writing this for readability):

[TABLE]

Next, we use

[TABLE]

After some algebra, we eventually get

[TABLE]

Thus, we obtain the formula we wanted

[TABLE]

The formula above has a very intuitive geometric meaning. The first term represents the first reflection from the top layer. The term $T_{+-}T_{-+}$ represents transmission into the next layer and back out. For two interfaces, $C$ will be $1$, but in general, it will be a reflection coefficients defined recursively using the next adjoining layer. The $R_{--}$ represents reflection from below the interface and the $C$ corresponds to a reflection from the next subsequent interface from above which in this case is at $r=R$. The exponent $p$ corresponds to how many such reflections occur before the wave transmits back to the upper layer. We will see later that terms such as $\frac{h_{-}^{(1)}}{h_{-}^{(2)}}$ correspond to travel times of each such interaction between two adjacent hypersurfaces (each being either a boundary or an interface).

Evaluating $1/D$

To proceed, we will use the asymptotic solutions to the ODE where the prefix $n$ means it satisfies the Neumann boundary conditions. Asymptotically, we must distinguish the various regimes for the different types of rays that may occur: reflecting, turning, grazing, gliding, evanescent, as well as combinations of these.

• •

Reflecting ($0<p<R/c(R)$): We use the linearly independent solutions to the ODE in the reflecting regime of the form (see [14, Appendix A])

[TABLE]

where $\delta_{\pm}$ is a function depending on $p$ that keeps track of phase changes for when a ray turns, and $R^{*}$ is the radius of the ray. For the reflecting regime where the ray never turns, then $\delta_{\pm}=0$ and $R^{*}=R$. For a general eigenfunction that does not necessarily satisfy the boundary conditions, we remove the subscript $n$ from the above definitions.

It is useful to note that in this reflecting regime, the transmission coefficients are independent of the frequency $\omega$ to leading order. Indeed, notice

[TABLE]

with similar formulas for the other terms $[2-],[1+],[1-]$. Then for $R_{++}$ (say), the $\omega$ in the first term above cancels from the numerator and denominator in the formula for $R_{++}$ so we have

[TABLE]

where $F$ is independent of $\omega$. Analogous results hold for the remaining reflection and transmission coefficients so we conclude

Lemma B.1**.**

To leading order as $\omega\to\infty$, the reflection and transmission coefficients are independent of $\omega$.

At the outer boundary $r=1$, when identifying $U_{n}$ with its principal term in the WKB expansion, the Neumann condition is $\partial_{r}U_{n}(1)=0$, which gives when $\omega\to\infty$

[TABLE]

Then

[TABLE]

Recall that for the calculation of $D$ that we need in (4.28), we replace the above $\omega_{n}$ by a general $\omega$ due to the contour integration and residue formula. Then for $\omega\to\infty$

[TABLE]

We finally obtain the expression for $r,r^{\prime}>b$

[TABLE]

It is convenient to set

[TABLE]

and let us label

[TABLE]

and

[TABLE]

where $R^{*}$ is the turning radius of the ray as in [14], and for the reflecting regime, $R^{*}=R$. Now $h^{(2)}_{+}(1)/h^{(1)}_{+}(1)=\exp(\operatorname{i}2\omega\Phi_{+})$. Observe that in the formula for $A$, $\frac{h^{(2)}_{+}}{h_{+}^{(1)}}=1$. We thus have $D^{-1}U_{l}(r)U_{l}(r^{\prime})$ is equal to

[TABLE]

where

[TABLE]

The first term in the sum has no $A$ in the coefficient since it represents going from $r$ to $r^{\prime}$ via $r=1$ with no interface interaction. The next two terms correspond to a direct ray from source and receiver located at radius $r$ and $r^{\prime}$. The fourth term represent a path from $r$ to $r^{\prime}$ via an interface interaction.

Next,

[TABLE]

To ease notation, note that terms like $\exp(-2\operatorname{i}\Phi_{\pm})$ correspond to 2-way radial travel times between two interfaces (one of which could be a boundary). So denote

[TABLE]

We then have

[TABLE]

Next observe that for a positive integer $q$

[TABLE]

Hence, the above sum becomes

[TABLE]

The boundary condition at $r=R$ forces $C=1$ here. We have shown that

[TABLE]

where $n_{(m_{0},m_{1},m_{2})}$ is a combinatorial coefficient.

Remark B.2**.**

We note that the combinatorial coefficients in the above formula have a special physical meaning. In [24], the author describes that in multilayered media, multiple waves travel different paths but arrive with the same travel times (kinematic analogs). Some of these can also have the same amplitude and phase characteristics (dynamic analogs). Hence, due to multiple scattering, there may be multiple waves with the same principal amplitude and travel time. If the corresponding ray is periodic, then all of these rays make a contribution to the trace and they are accounted for by the above combinatorial coefficient on the number of dynamic analogs for a particular ray. See [24, Figure 2] for examples of these dynamic analogs. The coefficients in the above formula agree with the simple counting argument in [24] for counting the number of dynamic analogs.

Radial travel times and amplitudes

Let us do the purely reflecting case first with no turning points since that is easier to index. Based on the above calculations, we want a convenient indexing to represent radial travel times and amplitudes of each wave constituent in the sum, and can be unified when we study the other regimes.

For $M=(m_{0},m_{1},m_{2})\in\mathbb{Z}_{\geq 0}^{3}$ let

[TABLE]

[TABLE]

Now we have corresponding amplitudes:

[TABLE]

Actually, we have to index more carefully since the amplitudes involving $A$ and $A^{2}$ above are not merely amplitudes but contain important phase information as well as seen in (B.11). Since this does not affect the main argument and only makes the index cumbersome, we opt not to do this.

Also, we note here that by forcing $C=1$, this enforces the inner Neumann boundary condition to leading order even on these generic solutions $U_{k}$. We choose to leave a generic $C$ in the formula since it will be needed for the case of multiple interfaces.

Substituting the Debye expansion for $\hat{G}$ in the 2

interfaces case

We will now insert the Debye expansion into the formula for $\hat{G}(r,r_{0},\omega)$ in (4.30).

First, we insert the leading order expansion (valid for $\text{Re}\ p>0$),

[TABLE]

to obtain, assuming $r,r_{0}>b$ (analogous formulas hold for $r>b,r_{0}<b$ or $r,r_{0}<b$ or $r<b,r_{0}>b$)

[TABLE]

The other regimes require only slight modifications to the computation above so we will be briefer on these since the notation can be unified to produce the above formula as well.

• •

Total internal reflection ($b/c_{-}(b)<p<b/c_{+}(b)$): In this case, we have a reflection from the interface with no transmission, which corresponds to an evanescent wave in $\Omega_{-}$.

Here $V_{-}$ will be evanescent with the form

[TABLE]

The reflection coefficients are computed identically except that in this case, $h_{-}^{(2)}=\exp\left(-\omega\int_{r}^{b}\left\lvert\beta_{-}\right\rvert\operatorname{d}\!r\right)$ so that the reflection and transmission coefficients are now complex valued.

The remaining formulas follow as before but are simpler since there is no propagation in the lower layer. Thus, we have the same formulas with $A$ replaced by $R_{++}$. With $r,r^{\prime}>b$, we have

[TABLE]

where

[TABLE]

and $\delta_{+}=0$ since there is no phase shift.

• •

Diving ($b/c_{+}(b)<p<1/c_{+}(1)$ or $R/c_{-}(R)<p<b/c_{-}(b)$): There are two possible cases, either the turning point is in the first layer (in this case, there will not be any reflections and transmissions so the analysis reduces to that of [14]) or the turning point is in the second layer, which requires further analysis than that of [14]. In the latter case, the rays in the first layer transmit into the second layer but then turn rather than reflect from $r=R$.

We summarize the WKB solution of (4.22) in the vicinity of a general turning point. A turning point, $r=R^{\star}$, is determined by

[TABLE]

Near a turning point, $r\approx R^{\star}$, and

[TABLE]

for an $q_{0}$ determined by a Taylor expansion. Away from a turning point,

[TABLE]

Matching asymptotic solutions yields

[TABLE]

From these one can obtain a uniform expansion, that is, the Langer approximation

[TABLE]

valid for $r\in[R,1]$. One obtains eigenfunctions corresponding with turning rays.

Up to leading order, where $r\gg R^{\star}$,

[TABLE]

B.2. Gliding and grazing cases

Recall that $u_{+}=S(A_{+,2,r}+AA_{+,1,r})$ If $u_{+}$ is an eigenfunction, then the surface boundary condition $u_{+}^{\prime}(1)=0$ determines the eigenvalues so

[TABLE]

and

[TABLE]

where $W(\cdot,\cdot)$ is the Wronskian. We then obtain

[TABLE]

Notice that $V_{-}$ is asymptotically [math] near $r=R$ and so the inner boundary condition is satisfied automatically.

With the above representation and following the ansantz in (B.1), one uses $C=1$, and after expressing the cosine term in $V_{-}$ in terms of complex exponentials, we use for $h^{(1)}_{-}$ and $h^{(2)}_{-}$

[TABLE]

Next, in this case there we have from (B.8) and (B.9) that $\delta_{+}=0$ (since the first layer has no turning points) while $\delta_{-}=-\pi/2$ due to the turning point phase shift as in [14]. The coefficent $A$ and $1/E$ are computed exactly as in the reflecting case except we have the extra $\pi/2$ phase shift coming from terms involving $V_{-}$. Hence, in the formula for $\hat{G}$, we have a contribution of $\pi/2$ to the KMAH index for each turning point in the ray path. More precisely, we will have $N_{M,i}=m_{1}+m_{2}+l_{i}$ where $l_{i}$ depends on the number of turning points. Hence, our earlier computation goes through where the KMAH index is the only difference.

One can see the phase shift analytically for each turning point along the ray in the formula for $A$. For example, a ray starting in the first layer that transmits into the second layer, if it transmits back to the first layer, the $\pi/2$ phase shift is accounted for in the terms $\tilde{T}_{-+}\tilde{T}_{+-}$ in (B.12). If it reflects from below the interface, the $\tilde{R}_{--}$ will have the $\pi/2$ shift.

• •

Reflection with gliding transmission ($p=b/c_{-}(b)$ and $b/c_{-}(b)<b/c_{+}(b)$ :

This is a case where a ray hits the interface at a critical angle so that there is a reflection but the transmitted ray begins tangent to the interface, and then propagates along the interface; inside $\Omega_{-}$, we are in the evanescent regime. In this case, $V_{+}$ is the same as before. Here we can also use

[TABLE]

while $V_{+}$ will be reflective

[TABLE]

from before.

Near $r=b$, we use the asymptotic formula $V_{-}\simeq 2T\pi^{1/2}q_{0}^{-1/6}\omega^{1/6}\mathrm{Ai}(-\omega^{2/3}(r-b))$ Let $S_{1}=2\pi^{1/2}q_{0}^{-1/6}\mathrm{Ai}(0)$ and $S_{2}=2\pi^{1/2}q_{0}^{1/6}\mathrm{Ai}^{\prime}(0)$. The transmission conditions (for $V$) then become with $d_{\pm}=\mu_{\pm}^{1/2}(b)$

[TABLE]

So then

[TABLE]

We have

[TABLE]

where $d=\omega^{5/6}d_{-}d_{+}^{-1}\beta_{+}^{-1/2}SS_{2}-\omega^{7/6}d_{+}d_{-}^{-1}\operatorname{i}\beta_{+}^{1/2}SS_{1}$. Hence, we have

[TABLE]

Similarly, $T\sim O(\omega^{-1/6})$ as $\omega\to\infty$.

In this case, to leading order we have

[TABLE]

Note that to leading order, this is analogous to the internally reflected case with no gliding. Hence, to leading order, only the travel time of the reflected ray can be detected and not the gliding portion.

• •

Grazing ($R^{*}=b\text{ or }R$): The detailed asymptotic analysis can be found in appendix C. For a grazing, turning ray, there are two new possibilities: The turning point is at the discontinuity $r=b$ or the turning point is at $r=R$.

Turning point at the discontinuity

Now we suppose $R^{*}=b$. Observe that due to the extended Herglotz condition, the lower layer becomes an evanescent regime where no propagation can occur, so using the asymptotic expansion of the Airy function, $V_{-}$ is exponentially decreasing and so $V_{-}=O(\omega^{-\infty})$ in $\Omega_{-}$. The inner boundary condition thus becomes automatically satisfied. Nevertheless, when restricted to the interface, $V_{-}\restriction_{r=b}$ is not [math] and is determined by the interface conditions. Hence, the interface conditions have the form

[TABLE]

for some $f_{1},f_{2}$ depending on $p$ and $\omega$. Next, note that $\mathrm{Bi}(0)=\sqrt{3}\mathrm{Ai}(0)$. Using (B.29) to represent $V_{+}$, the interface conditions have the form

[TABLE]

Away from the turning point, as in the turning point regime, $V_{+}$ is a linear combination of $\exp\left(\omega\int_{R^{\star}}^{r}\beta_{+}\,\operatorname{d}\!r^{\prime}-\pi/4\right)$ and $\exp\left(-\omega\int_{R^{\star}}^{r}\beta_{+}\,\operatorname{d}\!r^{\prime}+\pi/4\right)$ so upon replacing the factor $\sqrt{3}$ appearing in the equation for $B_{1}$ and $B_{2}$ with $e^{\operatorname{i}\pi/6}+e^{-\operatorname{i}\pi/6}$, this introduces an extra phase shift in the KMAH index, while the remaining calculations are analogous. Also, the formula for $A$ involving $h^{(2)}_{+}(1)/h^{(1)}_{+}(a)$ will also have this extra phase shift each time the ray turns due to the discontinuity.

In fact, we can calculate “reflection/transmission” coefficients (that is, an analog of them since there is no reflection/transmission in this case)too see what is happening. We use the ansatz

[TABLE]

for some $R$ to be determined. Since $\Omega_{-}$ is an evanescent regime, we could use the ansatz $V_{-}(r)=T\exp\left(-\omega\int_{r}^{b}\beta_{-}(r^{\prime})\operatorname{d}\!r^{\prime}\right)$. Then the same calculation for $R_{++},T_{+-}$ earlier would give to leading order as $\omega\to\infty$, $R=1$ and $T=0$ as expected since to principal order, there is no transmission to the other side.

We do the explicit computations since it is interesting to see how the discontinuity affects the leading order behavior. There are two calculations to do: Computing $D$ and computing $A$. We have

[TABLE]

To satisfy the Neumann condition, we evaluate at $r=1\gg b$ so we can use the leading order asymptotics of the Airy functions

[TABLE]

Writing the $\sin$ and $\cos$ terms in terms of complex exponentials and using the functions defined earlier $\Phi_{+}$ with $\delta=\pi/4$, we have

[TABLE]

The Neumann condition $\partial_{r}V_{n,+}(1)=0$ gives

[TABLE]

where we now use the actual eigenvalue $\omega_{n}$ and eigenfunction $V_{n,+}$. Thus, we have

[TABLE]

As before, we replace $\omega_{n}$ by the general $\omega$ and let $\omega\to\infty$

[TABLE]

Now we must compute $A$ which can be thought of now as the “reflection” coefficient from the interface. Near $r=b$, we can use the asymptotic formula

[TABLE]

It will be convenient to denote $S_{1}=2S\pi^{1/2}q_{0}^{-1/6}\mathrm{Ai}(0)$ and $S_{2}=-S2\pi^{1/2}q_{0}^{1/6}\mathrm{Ai}^{\prime}(0)\mu_{+,b}$. Below the interface, we use the ansantz $V_{-}=T\left\lvert\beta_{-}\right\rvert^{-1/2}\exp\left(-\omega\int_{r}^{b}\left\lvert\beta_{-}\right\rvert\operatorname{d}\!r\right)$ which satisfies the inner boundary condition to leading order.

Then the transmission conditions at $r=b$ are given by

[TABLE]

which leads to the matrix equation

[TABLE]

The determinant is $d=-\omega^{7/6}S_{1}\sqrt{3}\left\lvert\beta_{-}\right\rvert^{1/2}\mu_{-,b}-\omega^{5/6}S_{2}\left\lvert\beta_{-}\right\rvert^{-1/2}$ so the solution is

[TABLE]

So we obtain

[TABLE]

as $\omega\to\infty$, which modulo a normalization constant, the reflection coefficient is $1$ while $T=0$ to leading order, as expected.

We can now give a more explicit asymptotic formula for (B.28). First, note that $(-i+A)/(i+A)$ has modulus $1$ when using $A=-1/\sqrt{3}$ and angle $\arctan{-\sqrt{3}}=-\pi/3$. Hence, $(-i+A)/(i+A)=e^{-\pi/3}$, which to principal order, the extra phase shift the interface creates for the wave. Thus, we obtain

[TABLE]

Hence, in the earlier formula for $E$, we instead get

[TABLE]

where in the definition of $\Phi_{+}$, $\delta_{+}=\pi/6$, which is the adjusted phase shift from the turning ray that was $\pi/2$.

Turning point at $r=R$

It is possible that for certain turning rays, $R^{*}=R$, in which case the Neumann boundary condition $\partial_{r}V_{-}\restriction_{r=R=R^{*}}=0$ must be satisfied as well. This condition will be satisfied by using the representation (B.24) near the grazing point and also introducing $\mathrm{Bi}(x)$ in addition to $\mathrm{Ai}(x)$ above so that near $r=R^{*}$, $V_{-}$ is a linear combination of $\omega^{1/6}\mathrm{Ai}(-\omega^{2/3}q_{0}^{1/3}(r-R^{\star}))$ and $\omega^{1/6}\mathrm{Bi}(-\omega^{2/3}q_{0}^{1/3}(r-R^{\star}))$. So for $r\backsimeq R^{*}$, we use the ansantz

[TABLE]

Then

[TABLE]

Setting $\partial_{r}V(r)\restriction_{r=R}=0$ and using $-\mathrm{Ai}^{\prime}(0)=\mathrm{Bi}^{\prime}(0)/\sqrt{3}$ we get $C_{2}=\sqrt{3}$.

Now, $V_{+}$ is the same as the reflecting case, and the coefficient $A$ is computed identically. In the formula, the coefficient $C$ will change to account for the grazing at the boundary. To see this, for $r$ near $b$, we may write $V_{-}$ in terms of complex exponentials exactly as in the previous case where the turning point was at the interface. We obtain for $r$ near $b$

[TABLE]

Hence, in the formula for $A$ an $E$ in the previous sections, $C$ gets replaced by $e^{-\pi/3}$ which only affects the KMAH index each time the ray turns. Hence, the same computations as in the reflecting case go through with an adjusted KMAH index.

B.3. Multiple interfaces

We now consider an $N$-layered sphere whose wavespeed and density are smooth in each layer $j$ denoted $\Omega_{j}=\{d_{j}<r<d_{j+1}\}$. The index $j$ of the layers increases as $j$ increases. We have $N$ discontinuities $r=d_{j}$, $j=1,\dots,N$ and are also indexed in increasing radius. Hence $r=d_{N}$ is the surface and $r=d_{1}$ is the core-mantle boundary. As before, we consider the reflecting regime where there are no turning points. Analysis for the other regime will extend easily from the analysis we did in the two interface case.

• •

**Reflection ** ($0<p<R/c(R)$): We have the reflection coefficients $R_{j+1,jj+1}$ for a wave that reflects from interface $r=d_{j}$ from above. The first and last index indicate the wave began and ended in layer $\Omega_{j+1}$ while the middle index indicates which reflector the wave hit. Really, it is a function of $\omega$ and $l$. Then we have $T_{j+1,j}$ and $T_{j,j+1}$ as transmission from layer $j+1$ to $j$ and $j$ to $j+1$ respectively. Likewise, reflection from below interface $r=d_{j}$ is $R_{j,j,j}$. Corresponding to $C$ before, we will label $A_{j-1}$ corresponding to the total amplitude of outgoing waves at interface $r=d_{j-1}$.

Let $U_{j}=U\restriction_{\Omega_{j}}$. As before, we set

[TABLE]

[TABLE]

Then the same calculations as before lead to

[TABLE]

We denote

[TABLE]

and $M=(m_{1},\dots,m_{4(N-1)})$. We then define $Q_{M}:=Q^{M}$, that is, as the product of the amplitudes

[TABLE]

Note that $Q_{M}$ depends on $p$ but not $\omega$ using Lemma B.1. As before, we define $Q_{M,i}$ according to (B.17) with $A_{i}$ replacing $A$ in those formulas. The radial travel times $\tau_{M,i}$ will be constructed analogously using iteration from the two interface case. However, we do have to distinguish the different regimes to obtain the correct KMAH index.

For the radial travel times, define

[TABLE]

Here, $\delta_{j}$ depends on $p$, since if it is the reflecting regime for $\Omega_{j}$ (that is, the ray does not turn in $\Omega_{j}$, then $\delta_{j}=0$. If the ray turns in $\Omega_{j}$ but does not graze, then $\delta_{j}=\pi/2$. If it grazes, then $\delta_{j}=\pi/12$.

For $M=(m_{1},m_{2},\dots,m_{4(N-4)})\in\mathbb{Z}_{\geq 0}^{4(N-4)}$ let

[TABLE]

In such a case, iterating the calculation of (B.11) and (B.12) to obtain $E$ in (• ‣ B) from the single interface case, where $C$ gets replaced by $A_{j-1}$, we will have

[TABLE]

where $n_{M}$ is a combinatorial constant counting the number of dynamic analogs as in [24],

To simplify the notation, we assume that we are considering $U(r),U(r_{0})$ in layer $\Omega_{N}$ and $r_{N-1}=b$; the formulas are analogous for the other layers. As before we compute $D^{-1}U_{l}(r)U_{l}(r_{0})$ as

[TABLE]

As before, we denote the radial travel times as

[TABLE]

and corresponding amplitudes by

[TABLE]

As before, we would have to expand $A_{N}$ with a Neumann series, with each term contributing to the phase. However, the main form of the final formula does not change and so we opt not to do this in order to simplify the indexing.

• •

Turning/gliding/grazing/total internal reflection ($R/c(R)\leq p<1/c(1)$): There exists a minimal $l$ such that $d_{l}/c(d_{l})\leq p<d_{l+1}/c(d_{l+1})$. For the other regimes let $R^{*}$ be the turning radius of the deepest ray, and it depends on $p$. It is possible that the ray internally reflects at an interface or grazes it, in which case $R^{*}=r_{d_{l}}$. Due to Herglotz, all layers above $r_{d_{l+1}}$ are reflecting and we can repeat the analysis above to compute $U_{N},U_{N-1},\cdots,U_{l+3}$. Note that $\overline{\Omega}_{l+1}$ is where the ray turns, grazes, or internally reflects, so the region $r_{d_{j}}\leq r\leq r_{d_{j+2}}$ can be analyzed using the single interface case from before. To compute $U_{l+2},U_{l+1}$, we repeat the calculation of the single interfaces case to determine $A_{l+1}$ and $A_{l}$, where $A_{l}$ plays the role of $C$ in the two interface case and $A_{l+1}$ replaces $A$ in that case.

In all cases, we can unify the expressions so that to leading order we compute $\hat{G}(r,r_{0},\Theta,\omega)$, assuming $r,r_{0}>d_{N-1}$ (analogous formulas hold for the other intervals)

[TABLE]

It is important to note that $N_{M,i}$ depends on $p$, relating to a different phase shift from the various regimes described above.

B.3.1. Proof of Proposition 4.5: Wave trace near a gliding ray

Here, we will prove Proposition 4.5 showing the behavior of the wave trace near a gliding ray.

First, let $\gamma$ be a periodic ray with travel time $T$ that contains a gliding leg. We assume that other rays with travel time $T$ have the same number of reflected/transmitted legs or differ from $\gamma$ only through a rotation. Thus, there is an $\epsilon$ such that there are no periodic rays outside of $[\gamma]$ with travel time in $[T,\epsilon)$. We prove in section 2.2 that there is a sequence of nongliding, broken turning rays $\gamma_{m}$, $m=1,2,3,\dots$ converging to $\gamma$. Let $T_{m}$ being the travel time of these rays with ray parameter $p_{m}$. We would like to understand $\text{Tr}(\partial_{t}G)\restriction_{(T-\epsilon,T+\epsilon)}$.

Proof of Proposition 4.5.

First, let us assume there is only a single interface at $r=b$. When $\gamma$ hits the interface at a critical angle, the transmitted leg is tangent to the interface. As described in [8, p.181], when the angle of incidence is a little less than the critical angle, the ray of the transmitted wave has a turning point in the lower medium and later strikes the interface. It can be reflected from the interface (from below) and strike it again, and so on. Thus, the gliding wave, is a limit of waves which strike the interface from below $m=0,1,2,\dots$ times. These turning waves can be constructed with the standard WKB procedure we do in the turning regime. The limiting rays that strike the interface from below $m$ number of times are $\gamma_{m}$. There will be turning rays with travel times approaching $T$ from below that reflect from below the interface $m$ times. Following (4.35), the principal coefficient $a_{m}$ in the trace corresponding to this ray has the form

[TABLE]

where $C_{d}$ is independent of $m$, $Q_{m}$ is the product of the scattering coefficients, and the other quantities are explained there. Each term above remains bounded, but $Q_{m}$ and $\left\lvert\partial_{p}^{2}\tau_{m}\right\rvert^{-1/2}$ have decay properties that we will quantify as $m\to\infty$.

Since $\gamma_{m}$ enters the lower medium, reflects $m$ times, and exits into the upper medium, we have $Q_{m}(p_{m})=Q_{m}^{\prime}R_{m,--}^{m}T_{m,-+}$ for some uniformly bounded $Q_{m}^{\prime}$. We showed in Appendix B that for all non-gliding rays, the leading order contribution is $\int O(\omega^{3/2})e^{\operatorname{i}\omega(t-T_{m})}\operatorname{d}\!\omega$ while for $p=p_{G}$, the gliding ray, it is

[TABLE]

where $\epsilon>0$ is unknown, and it is even possible that $\epsilon=\infty$, which is essentially the case in [37, 15] albeit a slightly different setting. Hence, to leading order,

[TABLE]

We must make sure this sum is finite. First, we note in (B.5) that $R_{m,--}\to-1$ as $m\to\infty$. Next,

[TABLE]

Now, we already have $T_{m,+-}\to 0$ as $m\to\infty$ since $\beta_{m,-}\to 0$ but we need to know the rate this happens for the infinite sum above. Let $\Theta_{H}$ be the epicentral distance the gliding leg travels and $\Theta_{m,-}$ the epicentral distance of a turning segment. We know explicitly

[TABLE]

where $R^{*}_{m}<b$ is the turning radius. Next, we use that near the turning point, $r\approx R^{*}_{m}$, we have

[TABLE]

Hence,

[TABLE]

using that $R^{*}_{m}\to b$ as $m\to\infty$. We also have by our construction

[TABLE]

so for large $m$, $\beta_{m,-}=O(1/m)$ and hence $T_{m,+-}=O(1/m)$. Note that this is similar to estimate (6.17) in [8]. Also, the radial travel $\tau_{m}$ has the form $\tau_{m}=2\tau_{m}^{\prime}+2m\tau_{m,-}(b)$ where $\tau_{m}^{\prime}$ remains uniformly bounded. Hence, we obtain $\left\lvert\partial^{2}_{p}\tau_{m}\right\rvert^{-1/2}=O(1/\sqrt{m})$ (analogous to [5] and [8, Section 6.1]). Thus, the sum converges.

The same argument holds in the case of multiple interfaces. The limiting principal symbol $a_{m}$ will still involve a term of the form $T_{m,j,j-1}=O(1/m)$ where $r=d_{j}$ is the interface containing the gliding segment. In addition, the same argument above gives $\left\lvert\partial^{2}_{p}\tau_{m}\right\rvert^{-1/2}=O(1/\sqrt{m})$ which is all that is needed for a convergent sum. ∎

Appendix C Periodic grazing Ray

In this appendix, we will provide a more detailed analysis on the contribution of a periodic grazing ray to the trace formula. Our analysis closely follows [6, Chapter 1]. We do the analysis for $p$ near the grazing value $R/c(R)$ and then show the minor change necessary for $p$ near the value $b/c(b)$ corresponding to grazing at the interface. We will show that the leading order (as $\omega\to\infty$) contribution will have the “classic” form of (B.43) that can be handled with stationary phase while the lower order terms involve integrals of Airy functions where stationary phase does not apply. This is similar to the wave parametrix near a grazing ray described in [1] involving Airy functions.

We assume $U$ satisfies the inner boundary condition and $U_{n}$ satisfies both boundary conditions. We will need to compute

[TABLE]

We then replace $\omega_{n}$ by a general $\omega$. Using the asymptotic computation to sum the eigenfunctions computed earlier or using the computation in [40], we have the Green’s function representation

[TABLE]

Let $A_{r}$ and $B_{r}$ denote two linearly independent solutions to leading order for the equation (4.3) via solving (4.16) first. We will later pick

[TABLE]

and similarly for $B_{r}$ but using the Airy function $\mathrm{A}_{-}$, where $\mathrm{A}_{\pm}$ are Airy functions described in [5, 1]. Following similar notation as in section B and equation (B.1), we write $U_{n}$ restricted to the first layer $\Omega_{+}$

[TABLE]

for coefficients $S$ and $A$ that depend on $p$ and $\omega$, and $A$ was computed as (B.6). Similar to (B.1), for $U_{n}$ restricted to the second layer we set

[TABLE]

We do not add the $(\pm)$ superscripts for $A_{r},B_{r}$ since it will be clear in context based on which $r$ value we are evaluating. Note that $A$ is computed to be the same as (B.6) to satisfy the transmission conditions where $h_{+}^{(2)}=A_{b}$, $h_{+}^{(1)}=B_{b}$, and similarly for $h_{-}^{(1)},h_{-}^{(2)}$ in the formula. The Neumann inner boundary condition to leading order is

[TABLE]

so

[TABLE]

where a specific eigenvalue $\omega_{n}$ is being used, and for a radial function $D_{r}$, we use the notation $D^{\prime}_{b}=\frac{d}{dr}\restriction_{r=b}D_{r}.$ Thus, we get

[TABLE]

Since $U_{n}$ is an eigenfunction, then $\partial_{r}U_{n}=0$ at $r=1$ gives

[TABLE]

when $\omega=\omega_{n}$. Thus, we can write

[TABLE]

which implies

[TABLE]

where $W(A,B)$ is the Wronskian of $A_{r},B_{r}$ and is independent of $r$. We can now compute using (C.3)

[TABLE]

Thus,

[TABLE]

Note that even though $B_{s}^{\prime}=\frac{d}{dr}\restriction_{r=s}B_{r}$ and similarly for $A_{s}^{\prime}$, to leading order as $\omega\to\infty$, we have

[TABLE]

Next, following the computation in section B, the quantity $A^{k}$ above will consist of a sum of terms of the form

[TABLE]

where $A_{b^{\pm}}$ and $B_{b^{\pm}}$ comes from restricting $U^{(\pm)}$ to the interface $r=b$, and the last term comes from the quantity $C$ determined by the inner boundary condition. This last term is where stationary phase cannot be applied for $p$ near $R/c(R)$ while the other terms will be “classical” after using Airy function asymptotics. It will be convenient to use the multiindex $M=(m_{0},m_{1},m_{2},m_{3},m_{4},m_{5})\in\mathbb{Z}_{\geq 0}^{6}$.

When computing the trace $\int_{R}^{1}D^{-1}U(r)U(r)\rho dr$, we need to compute the quantities

[TABLE]

to leading order as $\omega\to\infty$. If these quantities are a symbol in $\omega$, as well as $B_{R}^{\prime}/A_{R}^{\prime}$ for $p$ near the grazing ray value $R/c(R)$, then we can just apply stationary phase to $(B_{1}^{\prime}/A_{1}^{\prime})^{k}$ using the asymptotic expansion of the Airy function as $\omega\to\infty$ by treating the rest of the integrand as the amplitude in the stationary phase calculation.

Thus, using (4.30) and the above computations, we get to leading order as $\omega\to\infty$

[TABLE]

where

[TABLE]

and

[TABLE]

is a symbol of order two, and $Q_{M}$ is a product of transmission and reflection coefficients described in appendix B. Let us write

[TABLE]

where

[TABLE]

We integrate by parts to obtain

[TABLE]

The first term is

[TABLE]

since $\mathrm{A}^{\prime}_{+}(\infty)/\mathrm{A}^{\prime}_{-}(\infty)=1$. This is the main term which has a classic form, where we can apply the method of steepest descent argument used in section 4.3. We just need to verify that the other term is indeed lower order.

After using the Airy equation, the second term becomes

[TABLE]

[TABLE]

where the subscript $R$ on $\chi_{R}$ means its evaluated at $r=R$ and

[TABLE]

Our integrand contains terms of the form

[TABLE]

so we use the substitution

[TABLE]

so $p=p(q)$ is a function of $q$ and we get

[TABLE]

Now we substitute

[TABLE]

to obtain

[TABLE]

Near the $p$ value $p_{g}:=R/c(R)$ corresponding to a periodic grazing ray is where stationary phase fails. If $p=p_{g}$, then $q=0$. Thus, we will do a Taylor series about $w=0$ and we have

[TABLE]

Applying the [6, proof of Proposition 9], the second term is indeed of order $\omega^{-2/3}$ and lower order than the principal term and can be disregarded. In fact, one can continue the Taylor expansion of the second terms and actually obtain lower order terms in the trace formula but we do not pursue this. Thus, taking the principal term gives us

[TABLE]

The analogous computation in [6, proof of Proposition 9], we have

[TABLE]

We are then left with

[TABLE]

For the other case where we consider periodic rays with a leg that grazes the interface, we need to do the above analysis for $p$ near $b/c(b)$. The above argument applies but the quantities $\left(\frac{A_{b^{+}}}{B_{b^{+}}}\right)^{m_{3}}$ and $\left(\frac{A_{R}^{\prime}}{B_{R}^{\prime}}\right)^{m_{5}}$ need to be interchanged in (C.5) and the rest of the argument below that.