The transmission problem in linear isotropic elasticity

TL;DR

This paper analyzes the behavior of elastic waves in a bounded domain with interfaces, revealing how wave transmission, reflection, and mode conversion depend on boundary conditions and geometric properties, and establishing unique determination of wave speeds.

Contribution

It provides a detailed microlocal analysis of elastic wave transmission across interfaces and derives conditions for unique wave speed recovery from boundary measurements.

Findings

Wave reflection, transmission, and mode conversion are characterized at interfaces.

Knott's equations are recovered in the context of elastic wave transmission.

Unique determination of P and S wave speeds from boundary data under convexity conditions.

Abstract

We study the isotropic elastic wave equation in a bounded domain with boundary with coefficients having jumps at a nested set of interfaces satisfying the natural transmission conditions there. We analyze in detail the microlocal behavior of such solution like reflection, transmission and mode conversion of S and P waves, evanescent modes, Rayleigh and Stoneley waves. In particular, we recover Knott's equations in this setting. We show that knowledge of the Dirichlet-to-Neumann map determines uniquely the speed of the P and the S waves if there is a strictly convex foliation with respect to them, under an additional condition of lack of full internal reflection of some of the waves.

Figures

Equations

$$[u]=0,[Nu]=0\text{on Γj, j=1,…,k},$$

$$c_{ijkl}=\lambda {\delta }_{ij}{\delta }_{kl}+\mu ({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}),$$

$$(Eu{)}_i={\rho }^{-1}\sum _{jkl}{\partial }_j{c}_{ijkl}{\partial }_l{u}_k,$$

$$u_{tt}-Eu=0,$$

$${\sigma }_{ij}(u)=\lambda \nabla \cdot u{\delta }_{ij}+\mu ({\partial }_j{u}_i+{\partial }_i{u}_j).$$

$$Nu:=\sum_{j}\sigma_{ij}(u)\nu^{j}\big{|}_{\partial\Omega},$$

$${\sigma }_p(-E)v=\frac{\lambda +\mu }{\rho }\xi (\xi \cdot v)+\frac{\mu }{\rho }|\xi {|}^{2}v,v\in {\mathbf{C}}^n,$$

$${\sigma }_p(-E)v=\frac{\lambda +2\mu }{\rho }\xi (\xi \cdot v)+\frac{\mu }{\rho }\left(|\xi {|}^2-\xi \xi \cdot \right)v.$$

$$c_p=\sqrt{(\lambda +2\mu )/\rho },c_s=\sqrt{\mu /\rho }$$

$$\{\begin{array}{ll}{u}_{tt}-Eu& =0\text{in R×Ω},\\ u{|}_{\mathbf{R}\times \partial \Omega }& =f,\\ u{|}_{t<0}& =0,\end{array}$$

$$(\Lambda f)_{i}=(Nu)_{i}=\sum_{j}\sigma_{ij}(u)\nu^{j}\big{|}_{\partial\Omega},$$

$$(Eu{)}_i={\rho }^{-1}({\partial }_i\lambda {\partial }_j{u}_j+{\partial }_j\mu {\partial }_j{u}_i+{\partial }_j\mu {\partial }_i{u}_j),$$

$$Eu={\rho }^{-1}(\mathrm{d}\lambda \delta u+2\delta \mu {\mathrm{d}}^{s}u),$$

$$({\mathrm{d}}^{s}u{)}_{ij}=\frac{1}{2}\left({\nabla }_i{u}_j+{\nabla }_j{u}_i\right),(\delta v{)}_i={\nabla }^j{v}_{ij},\delta u={\nabla }^i{u}_i,$$

$$\sigma (u)=\lambda (\delta u)g+2\mu {\mathrm{d}}^{s}u,$$

$$({\sigma }_p(-E)v{)}_i=\frac{\lambda +\mu }{\rho }{\xi }_i{\xi }^j{v}_j+\frac{\mu }{\rho }|\xi {|}_g^{2}v,v\in {\mathbf{C}}^n,$$

$$⟨\xi \times \eta ,\zeta ⟩=\omega (g^{-1}\xi ,g^{-1}\eta ,g^{-1}\zeta ),$$

$$g^{ij}{\alpha }_i{\zeta }_j=(\det g{)}^{-\frac{1}{2}}\det (\xi ,\eta ,\zeta ),$$

$$(\det g{)}^{1/2}{g}^{-1}\alpha =({\xi }_2{\eta }_3-{\xi }_3{\eta }_2,-{\xi }_1{\eta }_3+{\xi }_3{\eta }_1,{\xi }_1{\eta }_2-{\xi }_2{\eta }_1).$$

$$\xi \times \eta =(\det g{)}^{-\frac{1}{2}}g({\xi }_2{\eta }_3-{\xi }_3{\eta }_2,-{\xi }_1{\eta }_3+{\xi }_3{\eta }_1,{\xi }_1{\eta }_2-{\xi }_2{\eta }_1).$$

$$\nabla \times u=(\det g{)}^{-\frac{1}{2}}g({\partial }_2{u}_3-{\partial }_3{u}_2,-{\partial }_1{u}_3+{\partial }_3{u}_1,{\partial }_1{u}_2-{\partial }_2{u}_1).$$

$${\int }_D⟨Eu,v⟩\rho \mathrm{d}\mathrm{Vol}-{\int }_D⟨u,Ev⟩\rho \mathrm{d}\mathrm{Vol}={\int }_{\partial D}\left(⟨Nu,v⟩-⟨u,Nv⟩\right)\mathrm{d}A,$$

$$\int_{D}\langle Eu,v\rangle\rho\,\mathrm{d}\!\operatorname{Vol}=-\int_{D}\big{(}\lambda\langle\delta u,\delta v\rangle+2\mu\langle\mathrm{d}^{s}u,\mathrm{d}^{s}v\rangle\big{)}\mathrm{d}\!\operatorname{Vol}+\int_{{\partial}D}\sigma_{ij}(u)\nu^{j}v^{i}\,\mathrm{d}A,$$

$${\int }_{\partial D}⟨Nu,v⟩\mathrm{d}A.$$

$$(-Eu,u)={\int }_{\Omega }\left(\lambda |\delta u{|}^2+2\mu |{\mathrm{d}}^{s}u{|}^2\right)\mathrm{d}\mathrm{Vol},$$

$${\partial }_t(u_1,u_2)=\mathbf{E}(u_1,u_2):=(u_2,E{u}_1),(u_1,u_2){|}_{t=0}=(f_1,f_2).$$

$$\Vert (f_1,f_2){\Vert }_{\mathcal{H}}^2={\int }_{\Omega }\left(\lambda |\delta f_1{|}^2+2\mu |{\mathrm{d}}^s{f}_1{|}^2+|f_2{|}^2\right)\mathrm{d}\mathrm{Vol}.$$

$$(Nu{)}_j=\lambda (\delta u){\delta }_{j3}+\mu \left({\partial }_3{u}_j+{\partial }_j{u}_3-2{\Gamma }_{j3}^k{u}_k\right).$$

$$\begin{array}{rl}(Nu{)}_j& =\mu ({\partial }_3{u}_j+{\partial }_j{u}_3)-2\mu {\Gamma }_{j3}^{\nu }{u}_{\nu },j=1,2,\\ (Nu{)}_3& =\lambda ({\partial }_1{u}_1+{\partial }_2{u}_2)+(\lambda +2\mu ){\partial }_3{u}_3,\end{array}$$

$$({\partial }_t^2-c^2{\Delta }_g+A)u=0$$

Full text

First author partly supported by NSF Grant DMS-1600327

Second author partly supported by NSF Grant CMG-1025259 and DMS-1265958, and a Simons fellowship

Third author partly supported by DMS-1664683.

1. Introduction

The main goal of this work is to study the transmission problem in isotropic linear elasticity. Let $\Omega\subset\mathbf{R}^{3}$ be a smooth bounded domain. Let $\Gamma_{1},\dots,\Gamma_{k}$ be closed disjoint smooth surfaces (interfaces) splitting $\Omega$ into subdomains $\Omega_{k}$ with exterior boundary $\Gamma_{k-1}$ (with $\Gamma_{0}:=\partial\Omega$) and interior one $\Gamma_{k}$, see Figure 1, left. Assume that the density $\rho$ and the Lamé parameters $\mu$, $\nu$ are smooth up to those surfaces with possible jumps there. We also assume that at at every point, at least one coefficient has a non-zero jump. We impose the following transmission conditions

[TABLE]

where $[v]$ stands for the jump of $v$ from the exterior to the interior across any of those surfaces, and $Nf$ are the normal components of the stress tensor, see (2.3). We are motivated by the isotropic elastic model of the Earth where the density and the Lamé parameters jump across the boundary between the crust and the mantle, etc. We study the time-dependent elastic system, see (2.1).

The first goal of this paper is to describe qualitatively the microlocal behavior of solutions of this problem. At any interface $\Gamma_{i}$, an incoming S or P wave can generate two reflected waves, one S wave and one P wave through mode conversion and two transmitted ones. Then each branch can generate four more, etc., see Figure 1. In some cases, there might be a full internal reflection for one or both of the waves, and there could be no transmitted or reflected waves of a certain kind. In fact, the missing waves would be evanescent modes.

While works on geometric optics for the elasticity system exist (no transmission) [3, 13, 21, 22, 23, 31, 2], a comprehensive analysis of the transmission problem in linear elasticity has not been done to authors’ knowledge. In case of a flat surface and constant coefficients, some cases have been analyzed in the geophysics literature, see, e.g., [19, 20, 1, 25, 26]. In that case, if there is no full internal reflection, one looks for solutions in terms of potentials to reduce the number of variables; and the potentials of the four waves corresponding to an incoming one solve a system which decouples into a $4\times 4$ and a $2\times 2$ one, see also (7.5) and (7.6). Those equations were derived by Knott [16] and Zoeppritz [42] more than a century ago, see also [1]. In a recent paper [4], the hyperbolic-hyperbolic (HH) case is analyzed for variable $\lambda(x)$, $\mu(x)$ and $\rho=1$ but the construction for a curved boundary is partial only. The (HH) case is characterized by the wavefront of the Cauchy data on $\Gamma$: it could belong to projected S and P waves on either side of it, and in particular, there are no evanescent modes, see Section 5. This is just one of the many cases since we may have full internal reflection of some or both waves on one or both sides of $\Gamma$; and mode conversion to evanescent modes, see Section 9.3 for a summary. The most general study we are aware of is [40] where the coefficients are constant but cases other than the (HH) one are considered, even though not as extensively as we do it in this paper.

We analyze the general case of variable coefficients and a curved interface in all cases, away from glancing rays. We are interested in two main questions: is the problem well posed microlocally; and (control) can we create every configuration on one side with suitably chosen waves on the other. By doing that, we also compute the principal parts of the reflected and the transmitted waves. The microlocal well posedness reduces to showing the ellipticity of some $\Psi$DO system on $\Gamma$ with not particularly simple looking entries. Its solution serves as initial conditions for the corresponding transport equations for the hyperbolic of for the evanescent modes. In the flat, constant coefficients case, this system is actually the computation giving us the whole solution. Going back to the general case, in the (HH) microlocal region, we have four outgoing waves, each one being 3D vector-valued. This gives as a $12\times 12$ $\Psi$DO system for showing-well posedness. If we allow both S and P waves coming from both sides, we would have a $12\times 24$ system which we want to solve for some group of variables. The control question is reduced to solving the same system with a rearrangement of the unknowns: we are given the waves on one side and want to solve for the waves on the other.

Doing this analysis with brute force does not seem to be a promising approach. Instead, we look for inspiration in the geophysics (and the existing math) literature using the flat constant coefficient case as a starting point. We express the P and the S waves in terms of potentials, as the divergence and as the curl of such potentials on a principal symbol level first; and we extend this to an arbitrary order. We adapt this to the boundary value problem. Having such microlocal mode separation, we also split the S waves in the SV (shear-vertical) and SH (shear-horizontal) waves. This decomposition is valid on $\Gamma$ only, and depends on the point (and the codirection). Then we reduce those systems to more manageable decoupled $4\times 4$ plus $2\times 2$ ones for the outgoing solutions given the incoming ones; their extended versions are $4\times 8$ plus $2\times 4$ ones, see (9.12) and (9.13). If the boundary is flat and the coefficients are constant, those are exactly Knott’s equations [16]. Their ellipticity, needed to show well posedness, turns out to be a consequence of energy preservation (even though the determinant can be computed and analyzed [1]), another observation due to Knott. Ellipticity needed to show control can be verified easily and follows from the microlocal well posedness of the boundary value Cauchy problem.

We do this analysis in all microlocal cases with some or even all waves being evanescent; in that case we call them modes. The corresponding matrix symbols do not need to be recomputed; we just need to be careful which imaginary square roots to chose. Ellipticity based on energy preservation needs modifications though. Evanescent waves do not carry (high frequency) energy on the principal symbol level, at least.

We do such analysis for the boundary value problem for the outgoing solutions as well with Dirichlet or Neumann, homogeneous or not, boundary conditions. We also analyze the microlocal boundary value Cauchy problem. We start with the (principally scalar) acoustic equation first for two reasons: it is a needed ingredient in the analysis of the elastic system and SH waves behave as acoustic ones (no mode conversion).

We also study the surface waves propagating along the boundary (Rayleigh waves) or along an internal interface $\Gamma$ (Stoneley waves). Taylor [35] characterized Rayleigh waves as a propagation of singularities phenomenon when $n=2$ and $\partial\Omega$ is flat, and he also mentions that the analysis applies to the general case as well. The existence of such waves is due to lack of ellipticity of the Dirichlet-to-Neumann (DN) operator in the elliptic region on $\partial\Omega$ and in the elliptic-elliptic one on an internal interface. Restricted to the surface $\partial\Omega$ or $\Gamma$, they solve a real principal type of system; and the solution extends as an evanescent one in $\bar{\Omega}$. Yamamoto [40] viewed Stoneley waves in a similar fashion. A more detailed analysis of the Rayleigh and the Stoneley waves will appear in a work of Y. Zhang.

We also present an application of this analysis to the inverse problem of recovering the coefficients form the outgoing DN map. We recover first the lens relation associated with incoming S and P waves in the first layer $\Omega_{1}$; then we use the recent results by the authors [30] about local recovery of a sound speed (or a conformal factor) from localized travel times. By [2], we can recover $\rho$ in $\Omega_{1}$ as well, therefore we can recover all three coefficients $\mu$, $\lambda$ and $\rho$ there. In [30] we prove conditional Hölder stability as well which makes this approach for the inverse problem in this paper potentially stable as well; when it can be applied. In the case of no internal interfaces, this was done in [31]. The inverse problem for transversely anisotropic media is studied in [8]. The presence of interfaces however complicates the geometry considerably, see Figure 1 for the recovery of the coefficients in the deeper layers. The lens relation corresponding to a single S or P wave (ray) is multi-valued in general and there is no direct way to tell which branch is coming from which layer, roughly speaking. This makes the inverse problem much different. An essential difficulty following this approach is that there could be totally internally reflected rays in the interior side of one interface which never get out, not even through mode conversion. Then they cannot be generated by rays from the exterior (by “earthquakes”). We show that if there is no total internal reflection of S waves on the interface $\Gamma_{1}$ (from the interior), we can recover $c_{s}$ below it. This is more general than the result in [4] where $\rho=1$, and there is the implicit assumption that there is no full reflection of S and P waves. Since we do not recover all three coefficients below the first interface, we use arguments based on the geometry and the directions of the polarization only, which depend on the speeds only. Next, we also show that if there is no total internal reflection of $P$ waves as well, one can recover $c_{p}$ in $\Omega_{2}$. Those arguments can be used to get even deeper into $\Omega$ with the appropriate assumptions on the speeds.

2. Preliminaries

2.1. The elastic system

The isotropic elastic system in a smooth bounded domain $\Omega\subset\mathbf{R}^{3}$ is described as follows. The elasticity tensor is defined by

[TABLE]

where $\lambda$, $\mu>0$ are the Lamé parameters. Assume for now that the coefficients $\lambda$, $\mu$ and $\rho$ are smooth in $\bar{\Omega}$. The elastic wave operator is given by

[TABLE]

where $\rho>0$ is the density and the vector function $u$ is the displacement. The corresponding elastic wave equation is given by

[TABLE]

see, e.g., [26]. The stress tensor $\sigma_{ij}(u)$ is defined by

[TABLE]

Note that $Eu=\rho^{-1}\delta\sigma(u)$, where $\delta$ is the divergence of the 2-tensor $\sigma(u)$.

The Dirichlet boundary condition for $E$ is prescribing $u$ on the boundary; while the natural Neumann boundary condition is to prescribe the normal components of the stress tensor

[TABLE]

where $\nu$ is the outer unit normal on $\partial\Omega$. This is the operator appearing in the Green’s formula (2.15) for $E$ but also has the physical meaning as the infinitesimal deformation of the material in normal direction.

Let $\Gamma$ be a smooth surface where the coefficients $\rho$, $\lambda$, $\mu$ may jump. The physical transmission conditions across $\Gamma$ are the following. First, kinematic ones: the displacements $u$ on both sides of $\Gamma$ should match (no slipping of the material w.r.t. each other); and second, dynamical ones: the normal components $Nu$ on both sides should match (same traction). Therefore, if we declare one side of $\Gamma$ external and the other one internal, and denote by $[u]_{\Gamma}$ the jump of $u$ across $\Gamma$ from the exterior to the interior, we obtain the transmission conditions (1.1) on $\Gamma$. Note that in $[Nu]$, the operator $N$ depends on $\rho$, $\mu$ and $\lambda$ and has different coefficients on each side of $\Gamma_{j}$.

The operator $E$ is symmetric on $L^{2}(\Omega;\mathbf{C}^{3},\rho\,\mathrm{d}x)$. It has a principal symbol

[TABLE]

which can be also written as

[TABLE]

Taking $v=\xi$ and $v\perp\xi$, we recover the well known fact that $\sigma_{p}(-E)$ has eigenvalues $c_{p}^{2}$ and $c_{s}^{2}$ with

[TABLE]

of multiplicities $1$ and $2$ and eigenspaces $\mathbf{R}\xi$, and $\xi^{\perp}$, respectively. We have $c_{s}<c_{p}$. Those are known as the speeds of the P waves and the S waves, respectively. The eigenspaces correspond to the polarization of those waves. The characteristic variety $\det\sigma_{p}((\partial_{t}^{2}-E))=0$ is the union of $\Sigma_{p}:=\{\tau^{2}=c_{p}^{2}|\xi|^{2}\}$ and $\Sigma_{s}:=\{\tau^{2}=c_{s}^{2}|\xi|^{2}\}$, each one having two connected components (away from the zero section), determined by the sign of $\tau$.

Let $u$ solve the elastic wave equation

[TABLE]

with $f$ given so that $f=0$ for $t<0$ and all coefficients smooth in $\Omega$ (no transmission interfaces). The (outgoing) Dirichlet-to-Neumann $\Lambda$ map is defined by

[TABLE]

see (2.3), where $\nu$ is the outer unit normal on $\partial\Omega$, and $\sigma_{ij}(u)$ is the stress tensor (2.2).

2.2. An invariant metric based formulation

We have

[TABLE]

where we sum over repeating indices even if they are both lower or upper. This can also be written in the following divergence form

[TABLE]

where $\mathrm{d}^{s}u=(\partial_{j}u_{i}+\partial_{i}u_{j})/2$ is the symmetric differential, and $\delta=-(\mathrm{d}^{s})^{*}$ is the divergence of symmetric fields with the adjoint in $L^{2}$ sense.

To prepare ourselves for changes of variables needed in the analysis near surfaces that we will flatten out, we will write $E$ in an invariant way in the presence of a Riemannian metric $g$. We view $u$ as an one form (a covector field) and we define the symmetric differential $\mathrm{d}^{s}$ and the divergence $\delta$ by

[TABLE]

where $\nabla$ is the covariant differential, $\nabla^{j}=g^{ij}\nabla_{i}$, $u$ is a covector field, and $v$ is a symmetric covariant tensor field of order two. Note that $\mathrm{d}^{s}$ increases the order of the tensor by one while $\delta$ decreases it by one. Then we define $E$ by (2.10). We still have $\delta=-(\mathrm{d}^{s})^{*}$, where the adjoint is in the $L^{2}(\Omega,\mathrm{d}\operatorname{Vol})$ space of contravariant tensor fields, see, e.g., [24].

The stress tensor (2.2) is given by

[TABLE]

and then $Eu=\rho^{-1}\delta\sigma(u)$. The Neumann boundary condition $Nu$ at $\partial\Omega$ is still given by prescribing the values of $\sigma_{ij}(u)\nu^{j}$ on it as in (2.8). The operator $E$, defined originally on $C_{0}^{\infty}(\Omega)$ extends to a self-adjoint operator in $L^{2}(\Omega,\rho\,\mathrm{d}\!\operatorname{Vol})$. This extension is the one satisfying the zero Dirichlet boundary condition on $\mathbf{R}\times\partial\Omega$. In particular, this shows that the mixed problem (2.7) is solvable with regular enough data $f$ at least since one can always extend $f$ inside and reduce the problem to solving one with a zero boundary condition and a non-zero source term; and then use the Duhamel’s principle for the latter.

The principal symbol of $E$ in the metric setting is still given by (2.4) with the proper interpretation of the dot product there:

[TABLE]

where $\xi^{j}=g^{jk}\xi_{k}$ as usual. In particular, the speeds $c_{p}$ and $c_{s}$ remain as in (2.6). The eigenspaces of the symbol are still $\mathbf{R}\xi$ and $\xi^{\perp}$, the latter being the covectors normal to $\xi$. Notice that under coordinate changes, the coordinate expression for $u$ changes as well, as a covector.

We recall that the cross product on an oriented three dimensional Riemannian manifold is defined in the following way. If $\xi$ and $\eta$ are covectors at some fixed point $x$, then $\xi\times\eta$ is defined as the unique covector satisfying

[TABLE]

where $\langle\cdot,\cdot\rangle$ is the metric inner product of covectors, and $\omega$ is the volume form on the tangent bundle. To compute it in local coordinates, let $\alpha=\xi\times\eta$. Then we get

[TABLE]

where the latter is the determinant of the matrix with the indicated columns (also, the Euclidean volume form of them). Therefore, $(\det g)^{1/2}g^{-1}\alpha$ equals the Euclidean cross product

[TABLE]

This yields

[TABLE]

Similarly, the curl $\nabla\times u$ of a covector field $u$ is defined as the Hodge star of the exterior derivative $\mathrm{d}u$, and we have

[TABLE]

The divergence of $u$ is given by $\delta u=\nabla^{i}u_{i}$ and in particular, $\delta\nabla\times u=0$. We will use the notation $\nabla\cdot u$ for $\delta u$ as well.

One can verify that the double vector product of two covectors in the metric still satisfies $\xi\times\eta\times\zeta=\langle\xi,\zeta\rangle\eta-\langle\xi,\eta\rangle\zeta$, as in the Euclidean case.

2.3. Existence of dynamics

We assume now, as in the rest of the paper, that $\Omega$ can be expressed as a union of layers as explained in the Introduction and $\lambda$, $\mu$ and $\rho$ are smooth up to their boundaries with possible jumps at them. We also assume that $E$ is the metric based operator (2.10).

Lemma 2.1**.**

Let $\lambda,\mu,\rho$ be as above. Then $E$, defined originally on functions smooth up to $\Gamma_{1},\dots\Gamma_{k}$ and $\partial\Omega$, satisfying the transmission conditions (1.1), and zero boundary conditions on $\partial\Omega$, extends to a self-adjoint operator in $L^{2}(\Omega,\rho\,\mathrm{d}\!\operatorname{Vol})$.

Proof.

We start with Green’s formula. Let $D$ be a bounded domain with a smooth boundary so that $\lambda,\mu,\rho$ are smooth in $\bar{D}$. Then

[TABLE]

where $\mathrm{d}A$ is the area measure in ${\partial}D$ induced by $g$. To prove it, write

[TABLE]

since $Eu=\rho^{-1}\delta\sigma(u)$. The last integral equals

[TABLE]

Switch $u$ and $v$ and subtract the resulting formulas to prove (2.15).

Assume now that $u$ and $v$ are smooth up to the interfaces, may jump there and satisfy the transmission conditions (1.1). We apply (2.15) to $\Omega\setminus\Omega_{1}$, $\Omega_{1}\setminus\Omega_{2}$, …, $\Omega_{k}$ and sum up the results. Note that the outer normal to $\Omega\setminus\Omega_{1}$ at $\Gamma_{1}$ is the inner one at the same $\Gamma_{1}$ when viewed from $\Omega_{1}\setminus\Omega_{2}$, etc. As a result, we get (2.15) in $\Omega$ as well, despite the discontinuities because by the transmission conditions (1.1), all contributions from $\Gamma_{1},\dots,\Gamma_{k}$ cancel. By the zero boundary condition on $\partial\Omega$, the r.h.s. of (2.15) vanishes. Therefore, $E$ is symmetric.

To show that there is a natural self-adjoint extension, it is enough to show that the quadratic form $(-Eu,u)$ is bounded from below. For every smooth $u$ satisfying the Dirichlet boundary condition, by (2.10) we have

[TABLE]

which is non-negative.

We can write the Cauchy problem at $t=0$ for (2.1) with Dirichlet boundary conditions now as

[TABLE]

The operator $\mathbf{E}$ is self-adjoint on the energy $\mathcal{H}$ space with norm

[TABLE]

Then by Stone’s theorem, the Cauchy problem at $t=0$ for (2.1) with Dirichlet boundary conditions is solved by a unitary group. Problem (2.7) can be solved for regular enough $f$ by extending $f$ inside $\Omega$ and reducing it to a problem with a source but with homogeneous Dirichlet boundary conditions; and solving it by Duhamel’s formula. ∎

2.4. The Neumann boundary operator

Let $x=(x^{\prime},x^{3})$ be semigeodesic coordinates to a given surface $\Gamma$, with $x^{3}>0$ on one side of it, defining the orientation in the metric setup. The metric then takes the form $g$ in those coordinates with $g_{\alpha 3}=\delta_{\alpha 3}$ for $1\leq\alpha\leq 3$. Then, see also [31],

[TABLE]

Therefore,

[TABLE]

where $\nu=1,2$ and we used the fact that $\Gamma_{33}^{k}=\Gamma_{3k}^{3}=0$.

3. Geometric optics for the wave equation with $\Psi$DO  lower order terms

We recall the well known geometric optics construction for a hyperbolic pseudo-differential equation generalizing the acoustic wave equation, see, e.g., [36, 37]. We allow the equation to be a system but we still assume that the principal part is scalar, see also [9]. In this generality, the construction is done in [36, VIII.3]. We are not going to formulate results about the propagation of the polarization set which can be derived from [9]. The reason to do study the acoustic equation in this generality is two-fold. First, the elastic system decomposes into such pseudo-differential equations; and second, SH waves propagate like acoustic ones as we show below.

3.1. The Cauchy Problem with data at $t=0$

Our interest is in the acoustic wave equation with lower order classical pseudo-differential term $A\in\Psi^{1}$

[TABLE]

with Cauchy data $(u,\partial_{t}u)=(h_{1},h_{2})$ at $t=0$. Here, ${g}$ is a Riemannian metric that we include in order to have the flexibility to change coordinates easily; and $\Delta_{{g}}$ is the Laplace-Beltrami operator. The distribution $u$ is vector valued and $A$ is a matrix valued $\Psi$DO. Up to lower order terms, $c^{2}\Delta_{g}$ coincides with $\Delta_{c^{-2}g}$. The characteristic variety $\Sigma$ is given by $\tau^{2}=c^{2}|\xi|_{{g}}^{2}$ and has two connected components $\Sigma_{\pm}$ corresponding to $\tau<0$ and $\tau>0$, away from the zero section (notice the convention that $\tau<0$ corresponds to $\Sigma_{+}$). We are looking for solutions of the form

[TABLE]

modulo terms involving smoothing operators of $h_{1}$ and $h_{2}$, defined in some neighborhood of $t=0$, $x=x_{0}$ with some $x_{0}$. This parametrix differs from the actual solution by a smoothing operator applied to $\mathbf{h}=(h_{1},h_{2})$, as it follows from standard hyperbolic estimates. The signs $\sigma=\pm$ correspond to solutions with wave front sets in $\Sigma_{\mp}$, respectively as it can be seen by applying the stationary phase lemma.

Here, $a_{j,\sigma}$ are classical amplitudes of order zero depending smoothly on $t$ of the form

[TABLE]

where $a_{j,\sigma}^{(k)}$ is homogeneous in $\xi$ of degree $-k$ for large $|\xi|$. The phase functions $\phi_{\pm}$ are positively homogeneous of order $1$ in $\xi$ solving the eikonal equations

[TABLE]

Such solutions exist locally only, in general. While the principal symbol is the only one determining the eikonal equations and therefore the geometry, the subprincipal symbol in (3.1) depending on the principal one of $A$, affects the leading amplitude below.

Since the principal symbol of the hyperbolic operator in (3.1) allows the decomposition $-\tau^{2}+c^{2}|\xi|_{{g}}=(-\tau+c|\xi|_{{g}})(\tau+c|\xi|_{{g}})$, in a conic neighborhood of $\Sigma_{+}$, one can apply a parametrix of $D_{t}-c|D|_{{g}}$ to write (3.1) there as

[TABLE]

with $A_{+}$ of order zero and $u_{+}$ being the sum of the $\sigma=+$ terms in (3.2). This is the case studied in [36, VIII.3] with a more general elliptic $-\lambda(t,x,D)$ replacing $\mathrm{i}c|D|_{{g}}+A_{+}$, allowing $u_{+}$ to be a vector function, and $A_{+}$ to be matrix valued.

The main tool is the “fundamental lemma” allowing us to understand the action of a $\Psi$DO $P$ on $e^{\mathrm{i}\phi}a$ in terms of a homogeneous expansion in $\xi$, see [36, VIII.7] and [38]. The lemma remains true for principally scalar systems and it is used for such in [36].

We recall the construction of the amplitude. Let $u_{+}$ be as the first term in (3.2) with the indices there dropped, corresponding to $\sigma=+$. We seek the amplitude of the form $a=a_{0}+a_{1}+\dots$ as in (3.3) but the upper index $(k)$ is a lower one now.

The order two terms in the expansion of $(\partial_{t}-\mathrm{i}\lambda(t,x,D))u$ cancel because $\psi$ solves the eikonal equation (3.4) with the plus sign. Equate the order $1$ terms, we must solve

[TABLE]

where $\lambda=\lambda_{1}+\lambda_{0}+\dots$ is the expansion of $\lambda$ and they are evaluated at $\xi=\nabla_{x}\phi$. In our case, $\lambda_{1}=-c(x)|\xi|_{{g}}$, therefore, $\partial\lambda_{1}/\partial\xi=-cg^{-1}\xi/|\xi|_{{g}}$, which for $\xi=\nabla_{x}\phi$ yields $\partial\lambda_{1}/\partial\xi=-c{g}^{-1}\nabla_{x}\phi/|\nabla_{x}\phi|_{{g}}=c^{2}{g}^{-1}\nabla_{x}\phi/\phi_{t}$. Therefore, the vector field in (3.6) is proportional to the vector field $(\phi_{t},c^{2}{g}^{-1}\nabla_{x}\phi)$ which is the Hamiltonian covector field of the wave equation (3.1) on $\Sigma_{+}$ identified with a vector one, since the Laplacian there is the one associated with the metric $\tilde{g}:=c^{-2}{g}$. As it is well known, this is also the geodesic vector field of $\tilde{g}$ in the tangent bundle. The potential-like term in (3.6) involves $\lambda_{0}=-A_{+}$, see (3.5). Now, the transport equation (3.6) is a first order linear ODE along the bicharacteristics for the vector valued $a_{0}$ with a matrix valued zero order potential-like term. Given initial conditions at $t=0$, it is solvable as long as $\phi$ is well defined.

The higher order transport equations for $a_{1}$, $a_{2}$, etc., are derived in a similar way. They are non-homogeneous, with the same left-hand side but on the right we have functions computed in the previous steps.

We return to (3.2) now and look for $u$ as a sum of four terms as indicated here, each one of the type we described. We can use the Cauchy data to derive initial conditions for the transport equations, see e.g., [28], to complete the construction.

The integrals appearing in (3.2) are Fourier Integral Operators (FIOs) either with $t$ considered as a parameter, or as $t$ considered as one of the variables. In the former case, singularities of $(h_{1},h_{2})$ propagate along the zero bicharacteristics. More precisely, for every $t$,

[TABLE]

where $\mathbf{u}:=(u,u_{t})$, $\mathbf{h}=(h_{1},h_{2})$ and

[TABLE]

and for $(x,\eta)\in T^{*}\mathbf{R}^{3}\setminus 0$, $\gamma_{x,\eta}$ is the geodesic issued from $x$ in direction $\tilde{g}^{-1}\eta$.

On the other hand, considering $t$ as one of the variables,

[TABLE]

where

[TABLE]

In the analysis below, we will consider $C_{+}$ only.

The construction above can be done in some neighborhood of a fixed point $(0,x_{0})$ in general. To extend it globally, we can localize it first for $\mathbf{h}$ with $\operatorname{WF}(\mathbf{h})$ in a conic neighborhood of some fixed $(x_{0},\xi^{0})\in T^{*}\mathbf{R}^{3}\setminus 0$. Then $u$ will be well defined near the geodesic issued from that point but in some neighborhood of $(0,x_{0})$ in general. We can fix some $t=t_{1}$ at which $u$ is still defined, take the Cauchy data there and use it to construct a new solution. Then we get an FIO which is a composition of the two local FIOs each one associated with a canonical diffeomorphism, then so is the composition. Then we can use a partition of unity to conclude that while the representation (3.2) is local, the conclusions (3.7) and (3.8) are global. In fact, it is well known that both $\mathbf{h}\mapsto\mathbf{u}$ and $\mathbf{h}\mapsto\mathbf{u}(t,\cdot)$ with $t$ fixed are global FIOs associated with the canonical relations in (3.7) and (3.8).

In particular, if $\Gamma$ is a smooth hypersurface, and $\gamma_{x,\xi}(t)$ hits $\Gamma$ for the first time $t=t(x,\xi)$ transversely locally, then $\mathbf{h}\mapsto u|_{\Gamma}$ is an FIO again with a canonical relation as $C_{+}$ above but with $t=t(x,\xi)$ and $\dot{\gamma}$ replaced by its tangential projection $\eta^{\prime}:=\dot{\gamma}^{\prime}$. Notice that $\tau=-|\xi|_{\tilde{g}}<0$ for $C_{+}$ and $\tau=|\xi|_{\tilde{g}}>0$ for $C_{-}$. Also, $|\tau|<|\eta^{\prime}|_{\tilde{g}}$ with equality for tangent rays that we exclude; therefore, $\operatorname{WF}(u|_{\mathbf{R}\times\Gamma})$ is in the hyperbolic region, as defined below.

3.2. The boundary value problem for the acoustic equation

Let $\Gamma$ be a smooth hypersurface near a fixed point $x_{0}$ given locally by $x^{n}=0$. We take $x=(x^{\prime},x^{n})$ to be local semigeodesic coordinates. We define $\Omega_{\pm}=\{\pm x^{n}>0\}$ to be the “positive” and the “negative” sides of $\Gamma$. At the beginning, we work in $\Omega_{+}$ only and omit the superscript or the subscript $+$ from the corresponding quantities. For all possible solutions $u$ (not restricted to incoming or outgoing ones) with singularities not tangent to $\Gamma$, we want to understand how the Dirichlet data $f:=u|_{\mathbf{R}\times\Gamma}$ and the Neumann data $h:=\partial_{\nu}u|_{\mathbf{R}\times\Gamma}$ are related. Once we have this, we can understand microlocally the boundary value problems with either Dirichlet or Neumamn boundary conditions, or with Cauchy data.

The analysis depends on where the wave front set of the Cauchy data is. Let $(f,h)\in\mathcal{E}^{\prime}(\mathbf{R}\times\mathbf{R}^{n-1})$ be supported near some $(t_{0},x^{\prime})$. Then $T^{*}(\mathbf{R}\times\mathbf{R}^{n-1})\setminus 0$ has a natural decomposition into the hyperbolic region $c^{2}|\xi^{\prime}|_{{g}}<\tau^{2}$, the glancing one $\tau^{2}=c^{2}|\xi^{\prime}|_{{g}}$, and the elliptic one $c^{2}|\xi^{\prime}|_{{g}}>\tau^{2}$. Each one has two disconnected components corresponding to $\mp\tau>0$. We will recall the analysis in the $\tau<0$ component in more detail and will point out the needed changes when $\tau>0$. Also, we will not analyze (a neighborhood of) the glancing region; for that, see, e.g., [36] for a strictly convex boundary. We are looking for a parametrix of the outgoing solution $u$ of (3.1) with boundary data $f$, i.e., the solution with singularities propagating in the future only. Solutions with singularities propagating to the past only will be called incoming.

3.2.1. The outgoing and the incoming Neumann operators

If $u_{\textrm{out}}$ is the outgoing solution with boundary data $f$ with $\operatorname{WF}(f)$ in the hyperbolic region, we call the operator $\Lambda_{\textrm{out}}f=\partial_{\nu}u|_{\mathbf{R}\times\Gamma}$ the outgoing Neumann operator. Similarly we define the incoming Neumann operator by $\Lambda_{\textrm{in}}$. In those definitions, it is implicit that the solutions are defined in $\bar{\Omega}$ and $\nu$ is the unit normal exterior to it. i.e., $\partial_{\nu}=-\partial_{x^{n}}$. If we have $\Omega_{\pm}$ as above, we use the notation $\Lambda^{\pm}_{\textrm{in}}$, $\Lambda^{\pm}_{\textrm{out}}$ to denote the four Neumann operators with the convention that we preserve $\nu$ for $\Omega_{0}$, i.e., $\nu$ is interior for it. If the coefficients of the wave equation are smooth across $\Gamma$, we have $\Lambda_{\textrm{out}}^{+}=\Lambda_{\textrm{in}}^{-}$, $\Lambda_{\textrm{in}}^{+}=\Lambda_{\textrm{out}}^{-}$ up to smoothing operators. In the transmission problem below however, this is not the case.

3.2.2. Wave front set in the hyperbolic region $c^{2}|\xi^{\prime}|_{{g}}<\tau^{2}$

Assume that $\operatorname{WF}(f)$ is in the hyperbolic region with $\tau<0$. We are looking for a representation of $u$ of the form

[TABLE]

with a phase function $\phi$ and an amplitude $a$.

The phase function solves the eikonal equation in (3.4) with the plus sign but with a boundary condition on the timelike boundary $x^{n}=0$ now

[TABLE]

The choice of the positive square root reflects the assumption $\tau<0$. In the hyperbolic region, there are two solutions depending on the choice of the sign of $\partial_{x^{n}}\phi$ at $x^{n}=0$. It is easy to see that what corresponds to outgoing solutions is the positive choice

[TABLE]

We solve (3.10) with this condition locally. To construct the amplitude, we solve the same transport equations (3.6) as above but with initial condition $a=1$ for $x^{n}=0$, i.e., the principal part $a_{0}$ of $a$ is one there; and all others vanish.

The case $\tau>0$ is similar: we seek the solution in a similar way but the sign in (3.10) is negative. This does not change the construction.

Incoming solutions are constructed similarly. We choose the negative square root in (3.11). in particular we get that the outgoing and the incoming Neumann operators are $\Psi$DOs of order one with principal symbols equal to $\mp\mathrm{i}$ multiplied by (3.11), see also Proposition 3.1 below.

3.2.3. Wave front set in the elliptic region $c^{2}|\xi^{\prime}|>\tau^{2}$. Evanescent waves

We proceed formally in the same way but the problem here is that the eikonal equation has no real valued solution because the expression under the square root in (3.11) is negative. It may not even have a complex valued solution. This is a well known case of an evanescent mode described by a complex valued phase function (and amplitude). We follow [11], see also [36, VIII.4]. Since the construction in [11] is done for the Helmholtz equation with a large parameter and in [36, VIII.4] it is done for an elliptic boundary value problem, respectively, we need to do them in our hyperbolic case as well, even though the construction is essentially the same. We assume that $(t,x,\tau,\xi^{\prime})$ belong to a conically compact neighborhood, contained in the elliptic region, of a fixed point there. Plugging the ansatz in the elasticity equation, we use the “fundamental lemma” for complex phase functions in [38, X.4] to get an asymptotic expansion which formally look the same as in the hyperbolic case. We are looking for a solution of the eikonal equation (5.1) for $\phi$ up to an error $O(|x^{n}|^{\infty})$ at $x^{n}=0$ as a formal infinite expansion of the form

[TABLE]

where $\psi_{j}$ are symbols of order $1$. We denote this class by $\tilde{S}^{1}$, and by replacing the order $1$ by some $m$, we denote by $\tilde{S}^{m}$ the corresponding class. To avoid exponentially large modes, we require $\Im\phi\geq 0$. To construct the formal series, we first write the eikonal equation (3.11) in the form

[TABLE]

(note that there are no incoming/outgoing choices here) and then differentiate it w.r.t. $x^{n}$ at $x^{n}=0$. If such a solution exists, the error term would not affect those derivatives. We have

[TABLE]

To find the higher order derivatives, we write (3.12) in the form

[TABLE]

with $F(x,\eta)$ homogeneous in $\eta$ of order one. Then

[TABLE]

Since $\partial_{x^{n}}\phi$ is a symbol of order one, we prove the claim. Note also that $\Im\phi\geq x^{n}(|\tau|+|\xi|)/C$.

The next step is to solve the transport equations. Since they have complex coefficients, they may not be solvable exactly and we solve them up to an $O(|x^{n}|^{\infty})$ error as well. The rest is as in [11] and [36, VIII.4].

Proposition 3.1**.**

In the hyperbolic region, $\Lambda_{\text{\rm out}}$ and $\Lambda_{\text{\rm in}}$ are $\Psi$DOs of order one with principal symbols

[TABLE]

In the elliptic one, they are $\Psi$DOs of order one again with principal symbols

[TABLE]

We recall that $\partial_{\nu}=-\partial_{x^{n}}$ in the coordinates we used to compute the principal symbols. The expressions we got are invariant however. In both cases, the DN maps are elliptic. As shown in [36], they are elliptic even in the glancing region but they belong to a different class of $\Psi$DOs. The principal symbols of the Neumann operators on the negative side $\Omega_{-}$ are similar but with opposite signs.

3.2.4. The boundary value problem with Dirichlet data

The problem of constructing the outgoing solution $u_{\text{out}}$ with Dirichlet data on $\mathbf{R}\times\Gamma$ was solved above when $\operatorname{WF}(f)$ is either in the hyperbolic of the elliptic region. Similarly, we construct $u_{\text{in}}$. Notice that in the elliptic region, the construction is the same for both. In particular, we proved Proposition 3.1 by taking the normal derivatives of those solutions.

Next, we can construct a reflected wave. Assume we have an incoming solution $u_{\text{in}}$ with singularities hitting $\Gamma$ transversely. We want to construct a solution $u$ equal to $u_{\text{in}}$ for $t\ll 0$ satisfying $u=0$ on the boundary. Then $f:=u|_{\mathbf{R}\times\Gamma}$ has a wave front set in the hyperbolic region only. We construct the reflected wave $u_{R}$ as the outgoing solution with Dirichlet data $-f$. Then $u=u_{\text{in}}+u_{R}$ is the solution we seek.

3.2.5. The boundary value problem with Neumann data

Consider the outgoing solution $u_{\text{out}}$ with boundary data $\partial_{\nu}u=h$ on $\mathbf{R}\times\Gamma$. We reduce it to the Dirichlet problem above by inverting the DN map in $\Lambda_{\text{out}}f=h$. Since the latter is elliptic in the two regions we work in, this can be done microlocally. Then we solve a Dirichlet problem. We do the same for the incoming solution.

If we want to construct a reflected wave so that the solution $u$ satisfies $\partial_{\nu}u=0$, we need to solve $N_{\text{out}}f=-\partial_{\nu}u_{\text{in}}|_{\mathbf{R}\times\Gamma}$ which is possible since $N_{\text{out}}$ is elliptic. Having $f$, then we construct the outgoing solution with that Dirichlet data.

3.2.6. The boundary value problem with Cauchy data

We are looking for a microlocal solution $u$ of the acoustic equation (3.1) satisfying $u=f$ and $\partial_{\nu}u=h$ on $\mathbf{R}\times\Gamma$ with given $f$ and $h$ having wave front sets in the hyperbolic region first. The global Cauchy problem is over-determined because the singularities can hit the boundary again and therefore the Cauchy data have a structure (consisting of pairs in the graph of the lens relation); therefore prescribing them arbitrarily is not possible. On the other hand, one can construct a microlocal solution locally, when the wave front sets of $f$ and $h$ are localized in small conic sets excluding tangential directions, until the singularities hit the boundary again. We are looking for $u$ as a sum of two solutions $u=u_{\text{in}}+u_{\text{out}}$, one incoming and the other one outgoing. To determine the boundary values of the two solutions and to reduce the problem to section 3.2.4, we need to solve

[TABLE]

where $u_{\text{in}}$ and $u_{\text{out}}$ are the boundary values of those solutions.

Let $\operatorname{WF}(f,h)$ be in the hyperbolic region first. Then on principal symbol level, the leading amplitudes solve

[TABLE]

where $\xi_{3}$ is defined by (3.11). This in an elliptic system. This shows that the matrix valued operator in (3.16) is elliptic (if we reduce the order of the second equation to [math] by applying an elliptic $\Psi$DO of order $-1$). Therefore, the Cauchy data determine uniquely a decomposition into an incoming and an outgoing solution, locally. This reduces the problem to the one we solved in section 3.2.

If $\operatorname{WF}(f,h)$ is in the elliptic region, there is only one parametrix, no incoming or outgoing ones. The corresponding DN map $\Lambda$ is an elliptic $\Psi$DO of order one with principal symbol (3.15). Then for $(f,h)$ to be Cauchy data of an actual solution (up to smooth functions) it is needed that it belongs to the range of $(\textrm{\rm Id},\Lambda)$ (up to smooth functions). This makes this problem over-determined. If $h=\Lambda f$, a microlocal solution exists, as we showed above. It propagates no singularities away from $\Gamma$, and it does not propagate singularities along $\Gamma$ either (unlike the Rayleigh waves in elasticity which propagate along $\Gamma$).

3.3. The transmission problem

We recall the setup in section 3.2. We work locally in a small neighborhood of a point on $\Gamma$ and call one of its sides, $\Omega_{-}$ negative, the other one, $\Omega_{+}$, positive. For the speed $c$, we have $c=c_{-}$ in $\Omega_{-}$, and $c=c_{+}$ in $\Omega_{-}$, where $c_{-},c_{+}$ are smooth up to $\Gamma$ and $c_{-}\not=c_{+}$ pointwise. We impose the transmission conditions

[TABLE]

where $\nu$ is the normal derivative. Let $(x^{\prime},x^{3})$ be semi-geodesic coordinates near $\Gamma$ so that $\pm x^{3}>0$ in $\Omega_{\pm}$.

Let $u_{I}$ be an incident solution of the acoustic equation (3.1) with speed $c$ and background metric ${g}$ with a wave front set localized near a small conic neighborhood of some covector (at some time) approaching $\Gamma$ from the positive side. $\Omega_{+}$ As mentioned above, we consider singularities $(x,\xi)$ which move in the direction of $\xi$ only, i.e, associated with $\phi_{+}$ in (3.2), as we did in section 3. Then on $\operatorname{WF}(u_{I})$, with $t$ considered as a variable, we have $\tau<0$. Extend the speed $c$ form the negative to the positive side in a smooth way (recall that $c$ jumps across $\Gamma$) and extend $u_{I}$ smoothly across $\Gamma$ as a solution with that speed. Set

[TABLE]

Let $(x_{0},\xi_{0})$ with $x_{0}\in\Gamma$ be one of the singularities of $u_{I}$. We assume that $\xi_{0}$ is a unit covector w.r.t. $c_{+}^{-2}{g}$. We have that $\operatorname{WF}(f)$ is in the hyperbolic region $c_{+}|\xi^{\prime}|<-\tau$ in $\Omega_{+}$. We are looking for a parametrix $u$ near $x_{0}$ of the form

[TABLE]

where $u_{I}$ is incoming and restricted to $\bar{\Omega}_{+}$; $u_{R}$ is the reflected outgoing solution supported in $\bar{\Omega}_{+}$, and $u_{T}$ is the transmitted outgoing one or an evanescent mode, supported in $\bar{\Omega}_{-}$. It is enough to find the boundary values of those functions.

3.3.1. The hyperbolic-hyperbolic case

Assume that $\operatorname{WF}(f)$ is in the hyperbolic region in $\Omega_{-}$ as well, i.e., $c_{-}^{2}|\xi^{\prime}|^{2}<\tau^{2}$ on $\operatorname{WF}(f)$. If $c_{-}<c_{+}$ at $x_{0}$ (transmission from a fast to a slow region), that condition is satisfied regardless of $\xi_{0}^{\prime}$. If $c_{-}>c_{+}$(transmission from a slow to a fast region), existence of a transmitted ray depends on $\xi_{0}^{\prime}$. Let $\theta_{+}$ be the angle which an incoming ray makes with the normal, then the reflected angle will be the same and the angle $\theta_{-}$ of the transmitted ray, see Figure 2, is related to $\theta_{+}$ by Snell’s law

[TABLE]

which follows directly from (3.10) with $c=c_{-}$ and $c=c_{+}$ there, see also [29]. This relation shows that a transmitted ray will exist only if $\theta_{+}$ does not exceed the critical angle

[TABLE]

The transmission conditions (3.17) are equivalent to

[TABLE]

Assume now that we want to satisfy transmission conditions requiring continuity of $u$ and its normal derivative across the boundary. Then we get the following linear system for the leading terms $a_{T}^{(0)}$ and $a_{R}^{(0)}$ of the amplitudes $a_{T}$ and $a_{R}$:

[TABLE]

where

[TABLE]

In particular, this shows that the determinant of (3.23) is negative, and therefore, the system is solvable, i.e., elliptic after reducing the order of the second equation to zero. Since the system (3.22) is elliptic, it can be solved up to infinite order, i.e., we can find the all terms $a_{R,T}^{(k)}$ at $x^{n}=0$. The solutions serve as initial conditions for the transport equations of the corresponding modes.

Multiplying the first by the conjugate of the second equation, we get

[TABLE]

which can be considered (and justified) as preservation of the energy across $\Gamma$.

3.3.2. Total internal reflection

Assume now that $\operatorname{WF}(f)$ is in the elliptic region for $c_{-}$. This happens when $\theta_{+}>\theta_{\textrm{cr}}$. In that case, there will be no transmitted singularity. Indeed, we are looking for an evanescent mode in $\Omega_{-}$. Then $N_{\text{out}}^{-}$ in (3.22) is in the elliptic region. The analog of (3.23) then is

[TABLE]

where $\xi_{n}^{-}=\mathrm{i}|\xi_{n}^{-}|$ is pure imaginary and given by (3.15) times $\mathrm{i}$. Equivalently,

[TABLE]

Take the real part of the first equation multiplied by the conjugate of the second one to get

[TABLE]

In other words, on principal level, the whole energy is reflected and nothing is transmitted. We could have obtained this directly by solving (3.25), of course.

3.3.3. Incoming waves from both sides of $\Gamma$.

A more general setup is to assume incoming waves from each side, see Figure 2, right. We do not need to assume hyperbolic ones; they could be evanescent. In fact, this is an analogue of the Cauchy data case in the boundary value problem, see section 3.2.6. The point of view we adopt and will keep in the elastic case, is to classify the cases by the wave front set of the Cauchy data on the boundary.

We are interested in two questions: (i) well posedness of the transmission problem: given all incoming waves, is the problem well posed for the outgoing ones; and (ii) given all waves on one side of $\Gamma$, can we solve for all waves on the other one? We show that (i) is true as it can be expected (and well known). The answer to (ii) is not always affirmative; and when it is; this means that we can control the configuration on one side from the other one; in particular we can kill either the incoming or the outgoing wave on that side.

The hyperbolic-hyperbolic case. We assume now that the Cauchy data $(f,h)$ (the same on both sides by the transmission conditions) has a wave front set in the hyperbolic region on each side of $\Gamma$. Then on each side, we have two solutions: one incoming and one outgoing. Let $u^{+}_{\textrm{in}}$ and $u^{-}_{\textrm{in}}$ be the two incoming solutions from the positive and from the negative side, respectively, and let ${u}^{+}_{\textrm{out}}$, ${u}^{-}_{\textrm{out}}$ be the two outgoing ones. A usual, we assume no tangent rays. Then the transmission conditions are given by

[TABLE]

This is a generalization of (3.22) with one more wave added. If the corresponding principal amplitudes are $a^{+}_{\textrm{in}}$, $a^{-}_{\textrm{in}}$, $a^{+}_{\textrm{out}}$, $a^{-}_{\textrm{out}}$, we get

[TABLE]

Clearly, each matrix is elliptic. This implies that we have control from each side: given any choice of two amplitudes on one side, say $\Omega_{-}$, one gets an elliptic problem for finding the amplitudes on the other one, in this case $\Omega_{+}$.

We also get ellipticity for solving for the outgoing/incoming waves given the incoming/outgoing ones, i.e., the transmission problem is well posed. This also follows from energy conservation. Indeed, multiplying the first by the conjugate of the second equation, and then taking the real part above yields

[TABLE]

This energy preservation across the boundary implying in particular that if all incoming waves vanish, then so do the outgoing ones; i.e., that problem is elliptic.

The hyperbolic-elliptic case. We assume now that the Cauchy data $(f,h)$ (the same on both sides by the transmission conditions) has a wave front set in the hyperbolic region w.r.t. $c_{+}$ and in the elliptic one for $c_{-}$. Then in $\Omega_{+}$ we have two solutions: one incoming and one outgoing but in $\Omega_{-}$ there is only one (evanescent) solution. This case is analyzed in section 3.3.2 with $u_{\text{in}}^{+}=u_{I}$, $u_{\text{out}}^{+}=u_{R}$, $u^{-}$ (no incoming or outgoing ones) corresponding to $u_{T}$ there. We found out there that the incoming wave (or the outgoing one) determines uniquely the outgoing (respectively, the incoming) one and the evanescent one $u_{-}$. On the other hand, we cannot control $u_{\text{out}}^{+}$ and $u_{\text{in}}^{+}$ by choosing appropriately the evanescent mode $u^{-}=u_{T}$ appropriately; in fact $u_{\text{in}}^{+}$ alone determines the whole configuration already.

A slightly different point of view into this case is that we cannot have arbitrary (up to smooth functions) Cauchy data on $\Gamma$ in the hyperbolic region for $\Omega_{+}$, since that data falls in the elliptic region on the negative side, and then it has to be in the graph of the Neumann operator $\Lambda_{-}$. On other hand, if that data satisfy that compatibility condition, the solution in $\Omega_{+}$ consists of an incoming and a reflected wave. This is in contrast to the hyperbolic-hyperbolic case, where we can cancel one of the waves on the top, for example.

The elliptic-elliptic case. We assume now that the Cauchy data $(f,h)$ has a wave front set in the elliptic region w.r.t. both $c_{+}$ and $c_{-}$. It is interesting to see if we can have evanescent modes on both sides but still a non-trivial wave front set on $\Gamma$. We would need $(|\xi^{\prime}|_{g}^{2}-c_{+}^{-2}\tau^{2})^{1/2}=-(|\xi^{\prime}|_{g}^{2}-c_{-}^{-2}\tau^{2})^{1/2}$ which cannot happen. Therefore, there are no Rayleigh or Stoneley kind of waves in the acoustic case.

3.4. Justification of the parametrix

In each particular construction up to section 3.2.6, we constructed a parametrix satisfying the equation and the corresponding initial/boundary conditions up to a smooth error. Then the difference of the parametrix and the true solution satisfies all those conditions up to smooth errors. Standard hyperbolic estimates imply that the difference is smooth. In section 3.2.6, the Cauchy problem on a timelike boundary needs to be solved microlocally only and it is a tool to handle the transmission one. The justification of the parametrix for the latter can be done with the aid of [12, 39], guaranteeing smooth solutions if the transmission conditions (1.1) hold up to a smooth error only.

4. Geometric optics for the elastic wave equation

We study the Cauchy problem at $t=0$ and propagation of singularities in the elastic case. We present the geometric optics construction for the elastic wave equation in an open set first, where the coefficients are smooth. Such a construction is well known for systems with characteristics of constant multiplicities, see, e.g., [36, 37] and [9]. Our goal is to make the elastic case more explicit and to do a complete mode separation which we will use eventually near a boundary, see Proposition 4.1 below. The elastic case has been studied form microlocal point of view in [40, 21, 22, 23, 13, 3, 31].

Consider the elastic wave equation

[TABLE]

with Cauchy data $\mathbf{h}:=(h_{1},h_{2})$ at $t=0$. We want to solve it microlocally for $t$ in some interval and $x$ in an open set. The operator $E$ is associated with a Riemannian metric $g$ as in section 2.2. If $\lambda$, $\mu$ and $\rho$ are constant and $g$ Euclidean, one can use Fourier multipliers. In that case, let $\Pi_{p}=\Pi_{p}(D)$ be the projection to the p-modes, i.e., $\Pi_{p}$ is the Fourier multiplier $\hat{u}\mapsto(\xi/|\xi|)[(\xi/|\xi|)\cdot\hat{u}]$ and let $\Pi_{s}=\textrm{\rm Id}-\Pi_{p}$. It is easy to see that $\Pi_{s}$ is the Fourier multiplier $\hat{u}\mapsto-(\xi/|\xi|)\times(\xi/|\xi|)\ \times\hat{u}$. Also, we may regard $h=\Pi_{p}h+\Pi_{s}h$ as the potential/solenoidal (or the Hodge) decomposition of the 1-form $h$, see, e.g., [24]. Then , ${E}=c_{p}^{2}\Delta\Pi_{p}+c_{s}^{2}\Delta\Pi_{s}$. We have a complete decoupling of the system into P and S waves.

In the variable coefficient case, we will do this up to smoothing operators. We recall the construction in [36], which provides another proof of the propagation of singularities in this case. The principal symbol $\sigma_{p}(-E)$ of $-E$ has eigenvalues of constant multiplicities. Near every $(x_{0},\xi_{0})\in T^{*}\bar{\Omega}\setminus 0$, one can decouple the full symbol $\sigma(-E)$ fully up to symbols of order $-\infty$. Namely, there exist an elliptic matrix valued $\Psi$DO  $U$ of order [math] microlocally defined near $(x_{0},\xi_{0})$, so that

[TABLE]

modulo $S^{-\infty}$ near $(x_{0},\xi_{0})$, where the matrix is in block form; with an $1\times 1$ block on the lower right and a $2\times 2$ one on the upper left ($c_{s}^{2}\Delta_{g}+A_{s}$ is actually $c_{s}^{2}\Delta_{g}\,I_{2}+A_{s}$ with $I_{2}$ being the identity in two dimensions). Moreover, $A_{s}$ and $A_{p}$ are $\Psi$DOs of order one. In other words, the top non-zero block is scalar and the lower non-zero one is principally scalar. We recall this construction briefly. We seek $U$ as a classical $\Psi$DO with a principal symbol $U_{0}$ which diagonalizes $E$; there are many microlocal choices, and we fix one of them. Then

[TABLE]

where $R_{1}$ is of order one. Then we correct $U_{0}$ by replacing it with $U_{0}(I+K_{1})$ with some $\Psi$DO $K$ of order $-1$, i.e., we apply $I+K_{1}$ to the right and $(I+K_{1})^{-1}=I-K_{1}+\dots$ to the left to get

[TABLE]

where we used the fact that $(I-K_{1})R_{1}(I+K_{1})=R_{1}$ mod $\Psi^{0}$. Let us denote the matrix operator there by $G$. To kill the off diagonal terms on the right up to zeroth order, we need to do that for $GK-KG+R$. Note that $G$ and $K_{1}$ do not commute up to a lower order because they are matrix valued $\Psi$DOs. We look for $K$ in block form with zero diagonal entries and off-zero ones $K_{12}$ (an $1\times 2$ vector) and $K_{21}$ (a $2\times 1$ vector). If we represent $R_{1}$ in a block form as well, we reduce the problem to solving

[TABLE]

modulo $\Psi^{0}$. The solvability of this system on a principal symbol level follows by the general lemma in [36, IX.1] because $c_{s}\not=c_{p}$ but in this particular case, it is straightforward. Note that the principal symbols of $K_{12}$ and $K_{21}$ represent the coupling of the P and the S waves on a sub-principal symbol level, see also [3].

We apply $I-K_{2}$ to the left and $I+K_{2}$ to the right to kill the off diagonal terms of $(I+K_{1})^{-1}G(I+K_{1})$, etc. In fact, $U$ can be chosen to be unitary in microlocal sense [27]. In our case however, we prefer $U$ to be of order one.

From now on, we will do all principal symbol computation at a fixed point where $g$ is transformed to an Euclidean one (via the exponential map, for example) to simplify the notation. Then we will interpret the final result in invariant sense.

The principal symbol, of $U$, at that fixed point, will be chosen to be

[TABLE]

when $\xi_{3}\not=0$. The third column is the eigenvector $\xi$ associated with $c_{p}^{2}$, while the first and the second ones are a basis of the eigenspace of $\sigma_{p}(-E)$ associated with $\sigma_{s}$; and that basis is (micro) local only. In fact, a global one does not exist since those vectors are characterized as being conormal to $\xi$. In this particular case, we chose $\xi\times e_{1}$ and $\xi\times e_{2}$ with $e_{1}=(1,0,0)$, etc.

Recall that the principal symbol computations so far are at a single point where $g$ is Euclidean. To extend it to all points, an invariant way to choose $\sigma_{p}(U)$ is to replace the first and the second column there by $\xi\times e_{1}$ and $\xi\times e_{2}$ with $e_{1,2}$ considered as covectors, and the cross product as in (2.13). In other words, the first two columns in (4.5) are considered as vectors, then converted to covectors by the metric and multiplied by $(\det g)^{-1/2}$. Then we still get (4.9) but in $u^{s}$ we have curl in terms of the metric, see (2.14).

It then follows that microlocally, the elasticity system can be written as $(\partial_{t}^{2}-U^{-1}EU)w=0$ for

[TABLE]

where $w^{s}=(w_{1}^{s},w_{2}^{s})$ and $w^{p}$ is scalar. This system decouples into the wave equations

[TABLE]

with $A_{p,s}$ of order one, $R_{p,s}$ smoothing; the first one is a $2\times 2$ system and the second one is scalar. The first one has $\Sigma_{s}$ as a characteristic manifold, while the second one has $\Sigma_{p}$. Even though $U$ depends on the microlocal neighborhoods of the characteristic varieties $\Sigma_{s,p}$ we work in, the wave front sets of $U^{-1}f$, in those neighborhoods, we can apply the propagation of singularities results, or directly the microlocal geometric optics construction used below. Then we conclude that singularities in those neighborhoods propagate along the zero bicharacteristics of $\tau^{2}-c_{s}^{2}|\xi|^{2}$ and $\tau^{2}-c_{p}^{2}|\xi|^{2}$, respectively (which, of course, is well known). This implies a global result, as well.

For $u=Uw$ we get

[TABLE]

where $u^{s}$ and $u^{p}$ have wave front sets in $\Sigma_{s}$ and $\Sigma_{p}$, respectively. We call such solutions microlocal S and P waves. We have

[TABLE]

where $V_{p}$ and $V_{s}$ are of order zero and are formed by the lower order entries of $U$. Here $u^{s}$ can also be written as $u^{s}=D\times(w_{1}^{s},w_{2}^{s},0)+V_{s}w^{s}$.

Therefore, we proved the following.

Proposition 4.1** (mode separation).**

Let $u$ be a solution of the elastic wave equation in the metric setting in some open set in $\mathbf{R}\times\mathbf{R}^{3}$. Let $u^{p}$ and $u^{s}$ be $u$ microlocalized near $\Sigma_{p}$ and $\Sigma_{s}$, respectively. Then, microlocally, in any conic subset where $\xi_{3}\not=0$, there exist a scalar function $w^{p}$ and a vector valued function $w^{s}=(w_{1}^{s},w_{2}^{s})$ solving (4.7) so that $u=u_{p}+u_{s}$, where

[TABLE]

with $V_{p}$ and $V_{s}$ $\Psi$DOs of order zero and the curl in $D\times$ is in Riemannian sense.

The assumption $\xi_{3}\not=0$ does not restrict us. We can always rename the variables or rotate the coordinate system. On the other hand, the proposition does not provide a global mode separation. We are going to use it with $x^{3}$ being the distance to the boundary. Note also that $u$ and $w^{p}$, $w^{s}$ are related by (4.6).

In the geophysics literature, $w^{p}$ and $w^{s}$ such that $u^{s}=\nabla\times w^{s}$ (in our case, $w^{s}=(w_{1}^{s},w_{2}^{s},0)$) are called potentials. We have some freedom to choose $w^{s}$ so that (4.10) hold: adding an exact form to $(w_{1}^{s},w_{2}^{s},0)$ would not change the principal part of $u^{s}$ at least. One possible gauge to get unique $w^{s}$ is to take one of the components, in some coordinate system, to be zero. We have $w_{3}^{s}=0$ in (4.10). The analysis however must be restricted microlocally to $\xi_{3}\not=0$. In what follows, $x^{3}$ will be the normal coordinate to the boundary. Another choice is to require $w^{s}$ to be solenoidal, i.e., divergence free.

This proposition is a generalization of the well known representation of the solution of the isotropic constant coefficient elastic equation into potentials $u=\nabla w^{p}+\nabla\times w^{s}$ solving (4.7) with the operators $A_{p,s}$ and $R_{p,s}$ there vanishing. To guarantee uniqueness, it is often assumed that $w^{p}=(-\Delta)^{-1}\nabla\cdot u$, $w^{s}=-(-\Delta)^{-1}\nabla\times u$. We can prove a version of this in the variable coefficient case as well.

5. The boundary value problem for the elastic system. Dirichlet boundary conditions

Consider the elastic wave equation $u_{tt}-{E}u=0$ with boundary data $u=f$ on $\mathbf{R}\times\partial\Omega$. Assume that $f=0$ for $t\ll 0$ and we are looking for the outgoing solution, i.e., the one which vanishes for $t\ll 0$. We also introduce the notion of a microlocally outgoing solution along a single bicharacteristic requiring singularities of such a solution to propagate to the future. We define similarly incoming solutions by reversing time. Note that an outgoing solution does not need to consist of microlocally outgoing ones only since some incoming ones may be canceled at interfaces by outgoing ones. We will construct a parametrix of those solutions using the analysis in section (3.2). Moreover, we study the Cauchy data problem as well. We will use the analysis in the acoustic case essentially.

We work in semigeodesic coordinates $x=(x^{\prime},x^{3})$, with $x^{3}>0$ in $\Omega$. We denote the dual variables by $(\xi^{\prime},\xi_{3})$. The Euclidean metric then takes the form $g$ in those coordinates with $g_{\alpha 3}=\delta_{\alpha 3}$ for $1\leq\alpha\leq 3$. The analysis however works if we start with an arbitrary metric $g$ in $\mathbf{R}^{n}$, not just with the Euclidean one. Norms and inner products below are always in the metric $g$ or $g^{-1}$ (for covectors).

The phase space on the cylindrical boundary $\mathbf{R}\times\partial\Omega$ can be naturally split into the following regions (recall that $c_{s}<c_{p}$):

**Hyperbolic region: **

$c_{p}|\xi^{\prime}|_{g}<|\tau|$. Then $c_{s}|\xi^{\prime}|_{g}<|\tau|$ as well, so it is hyperbolic for both speeds.

**P-glancing region: **

$c_{p}|\xi^{\prime}|_{g}=|\tau|$. It is glancing for $c_{p}$ and hyperbolic for $c_{s}$.

**Mixed region: **

$c_{s}|\xi^{\prime}|_{g}<|\tau|<c_{p}|\xi^{\prime}|_{g}$. It is elliptic for $c_{p}$ but hyperbolic for $c_{s}$.

**S-glancing region: **

$c_{s}|\xi^{\prime}|_{g}=|\tau|$. It is glancing for $c_{s}$ and elliptic for $c_{p}$.

**Elliptic region: **

$|\tau|<c_{s}|\xi^{\prime}|_{g}$. Then $|\tau|<c_{p}|\xi^{\prime}|_{g}$, as well, so it is elliptic for both speeds.

We will not analyze wave fronts in the two glancing regions $|\tau|=c_{p}|\xi^{\prime}|_{g}$ and $|\tau|=c_{s}|\xi^{\prime}|_{g}$. For the purpose of the inverse problem, it is enough to analyze the propagation of singularities away from a set of measure zero. Therefore, there is no need to build a parametrix near the glancing regions (as in [32] or [41], for example) or work as in [12]; so we can avoid the glancing regions.

By the calculus of the wave front sets, the traces of microlocal P waves on $\mathbf{R}\times\partial\Omega$ have wave front sets in the hyperbolic region under the assumption that all singularities hit the boundary transversely. The traces of transversal microlocal S waves belong to $c_{s}|\xi^{\prime}|_{g}<|\tau|$, i.e, either to the hyperbolic, the mixed one, or to the p-glancing one. In particular, the trace of any solution of the elastic system with singularities hitting transversely, has wave front disjoint from the elliptic region. On the other hand, boundary values of solutions of the boundary value or the transmission problem may have wave front set on that surface, as Rayleigh and Stoneley waves do.

The analysis we have done so far, see next section, allows us to decouple the P and the S modes on the boundary completely by their polarizations. Then in terms of the potentials $w^{s}$ and $w^{p}$, we can think of the system as a decoupled one. When modes hit a free boundary, or a transparent one, however, the reflected and the transmitted modes may change type. The reason for this is that the boundary trace of an incoming S or P wave does not belong to the same subspace as that of an outgoing one.

5.1. Wave front set in the hyperbolic region

Let $f(t,x^{\prime})$ be supported near some $(t_{0},x_{0}^{\prime})\in\mathbf{R}\times\mathbf{R}^{2}$, where $\mathbf{R}^{2}$ represents $\partial\Omega$, flattened. Assume first that $\operatorname{WF}(f)$ is supported in the hyperbolic region. The later has two disconnected components determined by the sign of $\tau$ there. Let us assume that $\operatorname{WF}(f)$ is contained with the one with $\tau<0$; the $\tau>0$ case is similar. Then the characteristic varieties reduce to $\tau+c_{p}|\xi|_{g}=0$ and $\tau+c_{s}|\xi|_{g}=0$, respectively. We are looking for a parametrix of the outgoing solution of the form $u=Uw=u_{p}+u_{s}$ as in (4.8) with $w$ a potential. Note that this construction excludes $\xi_{3}=0$, which in our case corresponds to tangent rays which we avoid. We will work in a conic open microlocal region which does not contain such rays, i.e., $\xi_{3}\not=0$ there.

We seek the potentials $w^{p}$ and $w^{s}$ as geometric optics solutions as in section 3.2, i.e., of the form (3.9) (where the solution is called $u$, not $w$) with phases $\phi_{p}$ and $\psi_{s}$, respectively, and a scalar amplitude $a^{p}$ and a 2D vector-valued one $a_{s}=(a_{1}^{s},a_{2}^{s})$. The phase functions solve the eikonal equations

[TABLE]

and similarly for $\phi_{s}$, where $x^{\prime}=(x^{1},x^{2})$. The choice of the positive sign in front of the square root in the eikonal equation is determined by the choice $\tau<0$. By (4.10), the principal part of the amplitude of $u_{p}$ is $(D_{x}\phi_{p})a^{p}$ and that of $u_{s}$ is $D_{x}\phi_{s}\times(a^{s}_{1},a^{s}_{2},0)$. Restricted to the boundary, we have $\nabla_{x}\phi_{p}=(\xi^{\prime},\xi_{3}^{p})$, $\nabla_{x}\phi_{s}=(\xi^{\prime},\xi_{3}^{s})$, where

[TABLE]

We will use the notation

[TABLE]

Those are the codirections of the rays emitted from the boundary, see Figure 3. The angles $\theta^{p}$ and $\theta^{s}$ with the normal satisfy Snell’s law

[TABLE]

as it follows directly from (5.2), see also [29].

As we stated above, we are going to do all principal symbol calculations at $(t_{0},x_{0}^{\prime})$, where $g$ can always be arranged to be Euclidean.

In the hyperbolic region we work in, the expressions under the square roots are positive. The positive square roots guarantee that the singularities are outgoing. We determine next the boundary conditions for the transport equations. Since $u=Uw$, the boundary values of $w$ can be obtained from those of $u$ given by $f$ by an application of a certain $\Psi$DO. By the “fundamental lemma”, see [36, VIII.7] and [38], $Uw$ near the boundary is given by an oscillatory integral of the type (3.9) with the amplitude there multiplied by a classical symbol with principal part $U(x,\nabla_{x}\phi)$, where $\phi$ equals either $\phi_{s}$ or $\psi_{p}$ depending on which components of $w$ we take. Restricted to the boundary, we get

[TABLE]

with $U_{\textrm{out}}$ a classical $\Psi$DO  on $\mathbf{R}_{t}\times\mathbf{R}^{2}_{x^{\prime}}$ with principal symbol

[TABLE]

The subscript “out” is a reminder that we used the outgoing solution to define $U_{\textrm{out}}$. Similarly, we define $U_{\textrm{in}}$ using the incoming $u$. Its principal symbol is as above but with $\xi_{3}^{s}$ and $\xi_{3}^{p}$ having opposite signs. Note that $U$ acts locally in $\mathbf{R}_{t}\times\mathbf{R}_{x}^{3}$ while the two new operators act on $\mathbf{R}_{t}\times\mathbf{R}_{x}^{2}$. The symbol $\sigma_{p}(U_{\textrm{out}})$ is elliptic, in fact

[TABLE]

which also equals $\xi^{s}_{3}\langle\xi^{s},\xi^{p}\rangle$. The inverse of $\det\sigma_{p}(U_{\textrm{out}})$ is easy to compute and we do that below. To find the boundary conditions for $w=(w^{s}_{1},w^{s}_{2},w^{p})$, we write $w|_{x^{3}=0}=U_{\textrm{out}}^{-1}f$ (recall that all our inverses are parametrices). Then for $w^{p}$ and $w^{s}$ we get (4.9) with $\xi_{3}$ in all symbols replaced by $\xi_{3}^{p}$ for $u^{p}$ and $\xi_{3}^{s}$ for $u^{s}$. Once we have the boundary conditions for $w$, we construct $w$ near the boundary by the geometric optics construction (3.9). To get $u=u^{p}+u^{s}$, we apply $U$ to the result, see (4.8).

Remark 5.1**.**

In [23], Rachele showed that when $g$ is Euclidean, the leading amplitudes (polarizations) of $u^{p}$ and $u^{s}$ are independent of $\rho$ if we think of the three parameters being $(\rho,c_{s},c_{p})$ instead of $(\rho,\mu,\lambda)$. We will use this in Section 10.

In what follows, we will make the calculations above more geometric. By (4.10), $u^{s}$ and $u^{p}$ have representations of the kind (3.9) with the corresponding phase functions and matrix valued amplitudes having principal parts $f\to\xi\times(A^{s}f,0)$ and $f\to\xi A^{p}\cdot f$, where $A_{p}$ is a vector, and $A_{s}$ is a $2\times 3$ matrix. Then one can show that on the boundary, $h\mapsto\xi^{s}\times(A_{s}h,0)$ is the non-orthogonal projection to the plane $(\xi^{s})^{\perp}$ parallel to $\xi^{p}$, and $h\mapsto\xi^{p}A_{p}\cdot h$ is the non-orthogonal projection to $\xi^{p}$ parallel to the latter plane. In other words, they are the projection operators related to the direct sum $\xi^{p}\oplus(\xi^{s})^{\perp}$.

Finally in this section, we notice that the same analysis holds for the incoming solutions with given Dirichlet boundary data. Then in the formulas above, we have to take the negative square roots of $\xi_{3}^{p}$ and $\xi_{3}^{p}$ in (5.2).

5.2. Wave front set in the mixed region

Let $\operatorname{WF}(f)$ be in the mixed region next. We show below that the outgoing solution has a microlocal S wave only. The eikonal equation for $\phi_{s}$ still has the same real valued solution locally, corresponding to the outgoing choice of the solution $u_{s}$. On the other hand, the eikonal equation (5.1) for $\phi_{p}$ has no real solution. Indeed, we have $\nabla_{t,x^{\prime}}\phi_{p}=(\tau,\xi^{\prime})$ on $x^{3}=0$ and there is no real-valued function $\phi_{p}$ that could solve (5.1) and have such a gradient because in (5.2), $\xi_{3}^{p}$ would be pure imaginary. This is the case of an evanescent mode described in Section 3.2.3.

We are still looking for a solution of the form $u=u_{s}+u_{p}=U(w^{s},w^{p})$ but this time $w_{p}$, and therefore, $u_{p}$ is an evanescent mode as the one constructed in Section 3.2.3. The eikonal equation for $\phi_{p}$ implies, see (3.13), that $\xi_{3}^{p}$ in this case reduces to

[TABLE]

Then as in (5.5), (5.6), applying the “fundamental lemma” for FIOs with a complex phase, see [38, X.4], we deduce as before that the boundary values for $w$ are given by (5.5) with a classical $\Psi$DO $U$ having principal symbol as in (5.6) (with the new pure imaginary $\xi^{p}_{3}$). The operator $U$ is still elliptic because the determinant (5.7) has non-zero imaginary part. Then we can determine the boundary conditions for $w^{s}$ and $w^{p}$, construct the microlocal solutions, and apply $U$ to get $u$.

5.3. Wave front set in the elliptic region

Assume that $\operatorname{WF}(f)$ is in the elliptic region. Then we proceed as before, looking for both $w^{s}$ and $w^{p}$ as evanescent modes with complex phase functions. In this case, both $\xi_{3}^{p}$ and $\xi_{3}^{s}$ are pure imaginary with positive imaginary parts, see (5.8), and for $\xi_{3}^{s}$ we get

[TABLE]

We have

[TABLE]

Therefore, $U_{\textrm{out}}$ is elliptic and we can proceed as above and construct the solution as in Section 3.2.3.

5.4. Summary

We established that the Dirichlet problem is well posed microlocally and we have the following:

• (i)

$\operatorname{WF}(f)$ in the hyperbolic region: there are outgoing P and S waves.

• (ii)

$\operatorname{WF}(f)$ in the mixed region: there is an outgoing S wave only (plus an evanescent P mode).

• (iii)

$\operatorname{WF}(f)$ in the elliptic region: there are no outgoing waves; there are two evanescent modes.

6. The boundary value problem for the elastic system. Neumann boundary conditions and the Neumann operator

Assume now that we want to find the outgoing solution of the elastic wave equation with boundary data $Nu=h$. The strategy below is find the Dirichlet boundary data $f$ from this equation and then to proceed as in section 5. In other words, we want to solve $\Lambda f=h$ for $f$ microlocally if possible by showing that $\Lambda$ is elliptic (or not). Lack of ellipticity of $\Lambda$ in the elliptic region leads to Rayleigh waves, see, e.g., [5, 35, 33, 32].

6.1. Wave front set in the hyperbolic region

We are looking again for an outgoing solution of the type $u=u^{s}+u^{p}$ as in (4.8). The boundary values $w_{b}=w|_{x^{3}=0}$ of $w$ are computed by solving

[TABLE]

for $w_{b}$, compare with (5.5), where $\Lambda$ is the microlocalized Dirichlet-to-Neumann map (2.8), i.e., $\Lambda h:=Nu|_{x^{3}=0}$ for $u$ an outgoing microlocal solution of the elasticity equation with boundary data $u=h$ on $x^{3}=0$. We can use (6.1) and (5.6) to compute $\sigma_{p}(\Lambda)$.

We define the incoming $M_{\textrm{in}}$ in a similar way as in (6.1) but with $u$ being the incoming solution. More precisely, $M_{\textrm{in}}w_{b}$ is defined as $Nu|_{x^{3}=0}$ where $u$ is the incoming solution with boundary data $U_{\textrm{in}}w_{b}$. This also means that $M_{\textrm{in}}=\Lambda_{\text{in}}U_{\textrm{in}}$, where $\Lambda_{\text{in}}$ is defined as $Nu|_{x^{3}=0}$ with $u$ being the incoming solution. The operator $\Lambda$ the should be denoted by $\Lambda_{\text{out}}$ but we will keep the simpler one $\Lambda$. Below, we compute the principal symbols of $M_{\text{out}}$ and $M_{\text{in}}$. Combining that with (5.6), we can compute the principal symbol of $\Lambda$ as well but we will not need it.

By (2.16) and (5.6), in semigeodesic coordinates,

[TABLE]

Therefore,

[TABLE]

Similarly, we define $M_{\text{in}}$ to be the principal symbol of the same operator but related to the incoming DN map; i.e., the same as above but with $\xi_{3}^{s}$ and $\xi^{p}_{3}$ the negative square roots in (5.2).

A direct computation yields

[TABLE]

The determinant of $\sigma_{p}(M_{\text{in}})$ is the same. Since $U_{\textrm{out}}$ is elliptic, we get that $\Lambda$ is elliptic in the hyperbolic region as well. Therefore, we can invert $\Lambda$ microlocally and reduce the Neumann boundary value problem to the Dirichlet one, which can be solved as in section 5.1. More directly, we invert $\Lambda U_{\textrm{out}}$ and we get boundary conditions for $w$; which we use to solve the problem.

6.2. Wave front set in the elliptic region.

In this case, we seek both $w^{s}$ and $w^{p}$ as evanescent modes. The calculations are as in section 5 but $\xi_{3}^{s}$ and $\xi_{3}^{p}$ are pure imaginary as in (5.8) and (5.9). Then

[TABLE]

We have $c_{p}^{2}|\xi^{\prime}|^{2}-{\tau}^{2}>0$. For the third factor above, introduce the function

[TABLE]

Then, up to an elliptic factor, $\det\sigma_{p}(M_{\textrm{out}})$ equals $R(c_{s}^{-2}\tau^{2}|\xi^{\prime}|^{-2})$. It is well known and can be proven easily that on the interval $s\in(0,1)$, this function has a unique simple root $s_{0}$. This corresponds to $s_{0}c_{s}^{2}|\xi^{\prime}|^{2}=\tau^{2}$. Therefore, if we set $c_{R}(x)=c_{s}\sqrt{s_{0}}$, known as the Rayleigh speed, we get a characteristic variety

[TABLE]

on which (6.5) has a simple zero. Note that $0<c_{R}<c_{s}<c_{p}$. Since $U_{\rm out}$ is elliptic here, see Section 5.3, we get that $\Lambda$ is elliptic in the elliptic region away from $\Sigma_{R}$ and its principal symbol has a simple zero there. This generates the Rayleigh waves, see Section 8.2. For every $f$ with $\operatorname{WF}(f)$ in the elliptic region but away from $\Sigma_{R}$, we can proceed as above to solve the Neumann problem.

6.3. Wave front set in the mixed region

In this case, we seek both $w^{s}$ as a hyperbolic wave and $w^{p}$ as an evanescent one. The calculations are as in section 5 with $\xi_{3}^{s}$ real as in (5.2) and $\xi_{3}^{p}$ pure imaginary as in (5.8). Then $c_{p}^{2}|\xi^{\prime}|^{2}-{\tau}^{2}>0$ as well and for $\det\sigma_{p}(M_{\textrm{out}})$ we have an expression similar to (6.5) given, up to an elliptic factor, by $R(c_{s}^{-2}\tau^{2}|\xi^{\prime}|^{-2})$ with

[TABLE]

For $1<s<c_{p}^{2}c_{s}^{-2}$, which corresponds to the mixed region, $R$ is elliptic. This shows that, as above, one can construct $w|_{x^{3}=0}$ microlocally given $\Lambda f$. Then we construct $w^{s}$ and $w^{p}$, the latter as an evanescent mode; and then $u$. In particular, only microlocal S waves propagate from $\partial\Omega$.

6.4. Incoming solutions

The construction of incoming solutions (singularities propagating to the past only) is similar and we will skip the details. One can obtain them from the outgoing solutions by reversing the time.

7. The boundary value problem for the elastic system. Cauchy data

We analyze the boundary value problem for the elastic system on one side of $\Gamma$ with Cauchy data $u=f$, $\partial_{\nu}u=h$ on $\mathbf{R}_{t}\times\Gamma$. Similarly to section 3.2.6, we assume wave front set away from the glancing regions. This analysis is needed for the transmission problem when we want to control the behavior of the waves on one side by the other. We show in particular that this problem is well posed microlocally even though globally it is not, in general.

7.1. Wave front in the hyperbolic region

Assume first that the wave front set of $(f,h)$ is in the hyperbolic region. We are looking for a solution

[TABLE]

having both an incoming and an outgoing part, see Figure 4.

Then on $\Gamma$, we need to solve

[TABLE]

for the boundary traces $u_{\text{in},b}$ and $u_{\text{out},b}$ of $u_{\text{in}}$ and $u_{\text{out}}$. We pass to the corresponding solutions $w$ as in (6.1) to get

[TABLE]

Let $(a_{1,\text{in}}^{s},a_{2,\text{in}}^{s},a_{\text{in}}^{p})^{T}$ be the principal amplitude of $w_{\text{in}}$ and similarly for $w_{\text{out}}$. By the rotational invariance w.r.t. rotations in the $(\xi_{1},\xi_{2})$ plane (we justify this later), we can assume $\xi_{2}=0$. Then by (6.3),

[TABLE]

and similarly for $\sigma_{p}(M_{\text{in}})$, $\sigma_{p}(U_{\text{in}})$. Then on principal symbol level, (7.3) decouples into the following two systems

[TABLE]

and

[TABLE]

where

[TABLE]

We have

[TABLE]

This shows that the system (7.5) decouples to two $2\times 2$ systems after rewriting it as a system for the sum and the difference of the original vectors. The determinants of those two systems are $c_{s}^{-2}\tau^{2}\xi_{3}^{p}$ and $c_{s}^{-2}\tau^{2}\xi_{3}^{s}$, respectively; therefore, elliptic (after applying an elliptic operator of order $-1$ to the last two rows to equate their order with the rest, and we will use this notion of ellipticity below as well). Therefore, (7.5) is elliptic as well. Clearly, so is (7.6), which behaves as the acoustic case (3.16). Thus we proved the following.

Lemma 7.1**.**

The matrix valued symbol $(A_{\text{\rm in}},A_{\text{\rm out}})$ is elliptic.

Therefore, (7.3) is elliptic as well.

Lemma 7.1 remains true in the mixed and in the elliptic regions as well, where $\xi_{3}^{s}$ or $\xi_{s}^{p}$ could be pure imaginary as in (5.8), (5.9). Then there is no incoming/outgoing choice of the sign of $\xi_{3}^{s}$ and $\xi_{s}^{p}$ (which distinguishes $A_{\text{in}}$ and $A_{\text{out}}$) but this does not matter because later, we will multiply those expressions, when pure imaginary, with the “wrong” signs by zero, see (7.12), for example.

7.2. SV-SH decomposition of S waves

The principal amplitude of the S wave $u^{s}=D\times(w_{1}^{s},0,0)=(0,D_{3},-D_{2})w_{1}^{s}$ (plus smoother terms), see Proposition 4.1 and (4.9), corresponding to $w_{2}^{s}=0$, evaluated for $\xi_{2}=0$, has only its second component possibly non-zero. Then it is tangent to $\Gamma$ and normal to the direction of the propagation $\xi=(\xi_{1},0,\xi_{3})$ (as it should be because it is an S wave). In the geophysical literature (for constant coefficients and a flat boundary), such waves are called shear-horizontal (SH) waves since their polarization is tangent to the plane $\Gamma$. Equation (7.6) then describes the SH waves generated by the Cauchy data when $\xi_{2}=0$. Note that in our case, “horizontal” makes sense only at the boundary.

The $a_{2}^{s}$ terms appearing in (7.5) are the shear-vertical (SV) components of the potentials $w$ of the incoming and the outgoing waves. Indeed, using the subscript $b$ to indicate a boundary value (as we did above), when $w_{1,b}=0$, then the principal term of the outgoing/incoming $u^{s}_{b}$ is $(\mp\xi_{3}(x^{\prime},D^{\prime}),0,D_{1})w_{2,b}^{s}$, which gives us a principal amplitude perpendicular to the $\xi_{2}$ axis (and to the direction $\xi$ of propagation, of course). Then the oscillations happen in the $\xi_{1}\xi_{2}$ plane, vertical to $\Gamma$ (and parallel to $\xi$), hence the name. System (7.5) then describes how the SV and the P waves are created from given Cauchy data.

So far, the computations were done at a fixed point $x_{0}$ and a fixed covector $\xi^{0}$ at it, where the metric is chosen to be Euclidean. Then the orthogonal projection of the principal amplitude to $\Gamma=\{x^{3}=0\}$ (actually, to $T_{x_{0}}^{*}\Gamma$) is the SH component of it, while the projection to the plane through it and the normal is the SV component. We will do this decomposition microlocally near $(x_{0},\xi^{0})$ on the principal symbol level.

Note first that at $x_{0}$, there is a rotational invariance in the $\xi_{1}\xi_{2}$ plane. We already have a confirmation of that since we are free to choose coordinates in which $\xi_{2}=0$ and then we found out that the geometry of the rays and their principal amplitudes depend on the angles with the normal but not on $\xi$ in any other way. To derive this, we conjugate both symbols in (7.4) with the rotational matrix

[TABLE]

A direct computation yields

[TABLE]

at $x=x_{0}$. So far, we assumed that the metric was Euclidean at $x_{0}$. To get that, one can set $\tilde{\xi}=g^{-1/2}(x_{0})\xi$ which can be achieved by a linear change in the $x$ variables; then the Euclidean product in the $\tilde{\xi}$ variable corresponds to the metric one in the original $\xi$ one. Therefore, replacing $\xi$ above by $g^{1/2}(x_{0})\xi$ gives us the principal symbols in the original local coordinates. Varying the point $x_{0}$, we get principal symbols locally.

This allows us to define an SV-SH decomposition of S waves on a principal symbol level. In Proposition 4.1, if $u^{s}$ is the S wave of a solution with certain Cauchy data at $t=0$, then $u^{s}$ will be an SH wave on $\Gamma$ (up to lower order terms) if $\langle\nu,u^{s}\rangle|_{\Gamma}=0$ up to lower order terms applied to the Cauchy data, where $\nu$ is a unit normal covector field. It would be an SV wave on $\Gamma$ if $\langle\nu,(D\times u^{s})\rangle|_{\Gamma}=0$, up to lower order. An outgoing S wave $u_{\text{out}}^{s}$ near $\Gamma$, which is determined uniquely (up to a smooth term) by its Dirichlet data on $\Gamma$; and therefore by its potential $w_{\text{out},b}$ on $\Gamma$, is an SV wave on $\Gamma$, if $D^{\prime}\times w_{\text{out},b}=0$ up to a first order $\Psi$DO applied to $w_{\text{out},b}$, which corresponds to the requirement that the second component of $w_{\text{out},b}$ must vanish when $\xi^{\prime}=(\xi_{1},0)$. Here, $D^{\prime}$ is the tangential differential. To construct such SV waves, one can take the gradients on $\Gamma$ of scalar functions with non-trivial wave front sets. The $u^{s}$ wave is an SH one on $\Gamma$, if $D^{\prime}\cdot w_{\text{out},b}=0$ up to a lower order (divergence free). To construct such SH waves, one can take the curl on $\Gamma$ of scalar functions with non-trivial wave front sets.

7.3. Wave front in the mixed region

The P wave is evanescent, and there is only one (not incoming and an outgoing one). The number of the unknown amplitudes on the boundary is reduced by one, and the system can be seen to be over-determined. Indeed, then $\xi_{3}^{p}$ is pure imaginary and given by (5.8). We still define $A_{\textrm{in}}$ and $A_{\textrm{out}}$ as in (7.7). Then (7.5) becomes

[TABLE]

and (7.6) stays the same. By the expressions of the determinants following (7.8), the matrix $(A_{\text{in}},A_{\text{out}})$ is still elliptic in this case, i.e., Lemma 7.1 still holds. System (7.11) then is over-determined and solvable (uniquely) only if the r.h.s. belongs to a certain 3D subspace.

7.4. Wave front in the elliptic region

In this case, $\xi_{3}^{s}$ is pure imaginary as well as in (5.9), both waves are evanescent and the problem is overdetermined, as well. Equation (7.11) reduces to

[TABLE]

and Lemma 7.1 still holds with both $\xi_{3}^{p}$ and $\xi_{3}^{s}$ pure imaginary as in (5.8), (5.9); therefore we get an overdetermined system as well. In system (7.6), both amplitudes are equal and that system is overdetermined as well.

8. Reflection and mode conversion of S and P waves from a free boundary with Neumann boundary conditions

Let $\Gamma$ be a surface which separates an elastic medium from a free space (like the Earth from air). The natural boundary condition then is

[TABLE]

which means zero traction on $\Gamma$, i.e., no normal force, because the exterior has zero stiffness. We study reflection and mode conversion of S and P waves when they come from the elastic side of $\Gamma$ and hit $\Gamma$.

This is actually a partial case of the analysis of the boundary value problem with Cauchy data in Section 7 with zero Neumann and Dirichlet data. The strategy is the following. We take the trace $Nu_{I}$ of the incoming wave $u_{I}$ on the boundary and look for a reflected wave as a sum of an S and P wave as in (8.2) below. Then $Nu_{I}$ determines Neumann boundary conditions for those two waves. If $Nu_{I}$ has a wave front set in the hyperbolic region, we can recover the Dirichlet data for the reflected wave by inverting the elliptic $\Psi$DO $\Lambda U_{\textrm{out}}$ in (6.3). Knowing the Dirichlet data, we reduce the problem of constructing an outgoing solution as in section 5.1. If $\operatorname{WF}(Nu_{I})$ is in the mixed region, we use the construction in section 5.2. Finally, $\operatorname{WF}(Nu_{I})$ cannot be in the elliptic region since it corresponds to an incoming solution; therefore, Rayleigh waves cannot be generated by reflection of S and P waves. One can verify that the principal amplitudes of the reflected S and P waves can only vanish for a discrete number of incident angles (i.e., on a finite number of curves on the sphere of directions) because they depend analytically on $\xi$ and one can easily eliminate the scenario of one of the waves to vanish for all incoming directions. Those principal amplitudes can actually be computed and in the case of constant coefficients and a flat boundary, they have been computed in the geophysics literature, see, e.g., [1]. They do have zeros. For our purposes, it is enough to express their solution by Cramer’s Rule since we will prove that the determinant does not vanish. Vanishing amplitudes at finite number of angles is not an obstacle for the inverse problem we solve because the missing rays can be added to the data by continuity (but that may affect stability).

8.1. $\operatorname{WF}(u_{I,b})$ in the hyperbolic region

Assume that we have an incident P wave $u_{I}=u_{I}^{p}+u^{s}_{I}$, in other words a sum of microlocal solutions near $\Gamma$ with $\operatorname{WF}(u^{p}_{I})\subset\Sigma_{p}$ and $\operatorname{WF}(u^{s}_{I})\subset\Sigma_{s}$. As in Section 3, we will restrict the wave front set to $\tau<0$. We extend $u_{I}$ to a two sided neighborhood of $\Gamma$ as a microlocal solution by extending the coefficients $\lambda$, $\mu$ and $\lambda$ in a smooth way in the exterior. Set $u_{I,b}=u_{I}|_{\mathbf{R}\times S}$. It follows form the analysis above that $\operatorname{WF}(u_{I,b})$ is in the mixed region. As above, we assume no wave front set in the glancing region. In fact, $\operatorname{WF}(u_{I,b}^{p})$ is in the hyperbolic region while $\operatorname{WF}(u_{I,b}^{s})$ is there only if the angle of the corresponding rays with the normal is smaller than the critical one given by $c_{p}|\xi^{\prime}|=|\tau|$, and it is in the mixed one if the incident angle is greater than the critical one.

We look for a solution of the form

[TABLE]

where $u_{R}^{p}$ and $u_{R}^{s}$ are reflected P and S waves, respectively.

Let $x=(x^{\prime},x^{3})$ be semigeodesic coordinates near $x_{0}=0$ so that $x^{3}>0$ on the elastic side. All equalities below are at a fixed point $x_{0}$ which can be chosen to be [math] and modulo lower order terms for the amplitudes. As above, we assume without loss of generality that the metric $g$ is Euclidean at $x=0$ to simplify the notation. We can get the equations below by using (6.3). Let $w_{I}=(w_{1,I}^{s},w_{2,I}^{s},w_{I}^{p})$ and $w_{R}=(w_{1,R}^{s},w_{2,R}^{s},w_{R}^{p})$ be the solutions $w$ as in (4.8) related to $u_{I}$ and $u_{R}$. Since they solve (4.7), each singularity of the S or the P part of $w_{I}$ reflects by the laws of geometric optics. On the other hand, if $\theta^{p}$ is the angle which an incoming P singularity makes with the normal, then the corresponding angle $\theta^{s}$ of the reflected S singularity, see Figure 5, is related to $\theta^{p}$ by Snell’s law (5.4) as it follows directly from (5.2), see also [29]. Also, the incoming and the outgoing directions, and the normal belongs to the same plane, which determines the reflected direction uniquely. The same law applies to an incoming S wave generating a reflected P one. In the latter case, there is a critical incoming angle $\theta_{\text{cr}}=\arcsin(c_{s}/c_{p})$ of an S wave so that if $\theta^{s}>\theta_{\text{cr}}$, (5.4) has no solution for $\theta^{p}$. Then a reflected P wave does not exist and instead we have an evanescent mode, as we show below.

We need to solve

[TABLE]

for $w_{R,b}$. Since $M_{\text{out}}$ is elliptic in the hyperbolic region, (8.3) is microlocally solvable. We only need to verify that $w_{R}$ has non-trivial S and P components for almost all incoming rays.

We express $w_{I}$, $w_{R}$ and $w_{T}$ in the form (3.9) with phase functions solving (3.10) with for either $c_{p}$ or $c_{s}$ and a choice of the square root sign corresponding to the incoming or the outgoing property of each wave. The corresponding principal amplitudes are $(a_{1}^{s},a_{2}^{s},a^{p})$ subindices $I$, $R$, and $T$ distinguishing between the three waves.

Without loss of generality, we may assume $\xi_{2}=0$ as in section 7. We get, see (7.4),

[TABLE]

The system (8.4) is uniquely solvable, as we know. We determine $a_{1,R}^{s}=-a_{1,I}^{s}$ first, which says that the SH wave $U(a_{1,R}^{s},0)$ just flip a sign at reflection. The first two equations can be solved to get $a_{2,R}^{s}$ and $a^{p}_{R}$. If $a_{R}^{p}=a_{1,R}^{s}=0$, then $U(0,a_{2,R}^{s})$ is the SV wave oscillating in the plane normal to the boundary.

Let $w_{I}$ be a purely P wave, i.e., $w_{1,I}^{s}=w_{2,I}^{s}=0$. We want to find out when there is no reflected either P or an S wave. One could just solve the system but we will analyze it without solving it. If there is no reflected P wave, i.e., if $w^{p}_{R}=0$, then (8.4) implies that both components of the reflected wave must vanish as well which is a contradiction, unless $2\mu\xi_{1}^{2}-\rho\tau^{2}=0$, i.e., if $2c_{s}^{2}\xi_{1}^{2}=\tau^{2}$. This may or may not be in the hyperbolic region and defines a cone of incoming directions when it does. Now, assume that there is no reflected S wave, i.e., $w_{1,R}^{s}=w_{2,R}^{s}=0$. This is possible only when $\xi_{1}=0$, i.e., when the incoming P wave is normal to the boundary.

Now, assume that $w_{I}$ is an S wave. If there is a reflected S wave only, we are in the situation above with the time reversed — it can only happen for normal rays. Similarly, if there is a reflected P wave only, this can only happen for incident directions on a specific cone, or it does not.

8.2. Wave front set in the elliptic region, Rayleigh waves

We are looking for microlocal solutions satisfying $Nu=0$ with wave front set on the boundary in the elliptic region. We follow Taylor [35], where the coefficients are constant and $n=2$ but as noted there, the construction extends to the general case; and will sketch that extension. As shown in Section 6.2, $\Lambda$ has a characteristic variety $\Sigma_{R}$, see (6.6) and the determinant of its principal symbol, up to an elliptic factor near $\Sigma_{R}$, is given by $H:=\tau^{2}-c_{R}^{2}|\xi^{\prime}|^{2}$. Therefore, microlocal solutions to $Nu=0$ with boundary wave front sets on $\Sigma_{R}$ would solve a $\Psi$DO system on $\mathbf{R}_{t}\times\Gamma$ of real principal type in the sense of [9]. Here, $|\xi^{\prime}|$ is the norm of the covector $\xi^{\prime}$ in the metric on $\Gamma$ induced by $g$ (the latter is Euclidean in the isotropic elastic case). One can impose Cauchy data at $t=0$ to get unique (in microlocal sense) solution. Singularities propagate along the null bicharacteristics of $H$, i.e., along the null bicharacteristics of a wave equation on $\mathbf{R}_{t}\times\Gamma$ with speed $c_{R}$.

Next, one uses the solution on $\mathbf{R}_{t}\times\Gamma$ constructed above as Dirichlet data for a solution near $\Gamma$, in $\Omega$, as in Section 3.2.3.

8.3. Wave front set in the mixed region

This can only happen if there is a non-zero incident S wave hitting the boundary at an angle (with the normal) greater than the critical one $\theta_{\textrm{cr}}$, see (3.21). We are still looking for a solution of the kind (8.2), where $u_{I}^{p}=0$ and all singularities of $u_{I}^{s}$ hit the boundary at angles greater than $\theta_{\textrm{cr}}$. Then $u_{R}^{p}$ would be actually an evanescent mode (not actually a P wave by our definition because it would be smooth away from $\Gamma$). To find the boundary values for $w_{R}$, we need to solve (8.3) again with $M_{\text{out}}$ as in (6.3) but $\xi_{3}^{p}$ is given by (5.8). The matrix $M_{\text{out}}$ is elliptic, see (6.7). Once we have the boundary values for $w_{R}$, we can construct both solutions as in section 5. We also see that the reflected S wave cannot have zero amplitude except for possibly one incident angle; the proof is like in the hyperbolic case.

8.4. Summary

• (i)

An incident P wave produces a reflected P wave and a reflected S wave.

• (ii)

An incident S wave produces a reflected S wave. It produces a reflected P wave only if the incident angle is greater than the critical one; otherwise there is an evanescent P solution.

• (iii)

By time reversal, given an outgoing P wave, there are incoming S and P ones which produce that P wave and no S wave. The roles of those waves can be reversed only when the incident angle of the S wave is greater than the critical one.

• (iv)

An incident SH wave produces a reflected SH wave only.

9. The transmission problem for the elastic system

9.1. Transmission and reflections of incoming S and P waves. Zoeppritz’ and Knott’s equations

We are interested first how an incoming wave, either and S or an P one, is reflected and transmitted across $\Gamma$. We assume first that the wave front set of the incoming waves on the boundary is in the hyperbolic region on the other side of $\Gamma$ as well. This is a classical case with a long history. As in section 3.3, we assume that $\Gamma$ divides $\mathbf{R}^{3}$ locally into $\Omega_{+}$, where the waves come from, and $\Omega_{-}$, where they may transmit. Let, as above, $u_{I}$ be a microlocal solution of the elastic system. Similarly to (3.19), we are looking for a local solution of the form

[TABLE]

where the expressions in each parentheses is a decomposition into P and S waves, $u_{T}$, is supported in $\bar{\Omega}_{+}$, and $u_{I},u_{R}$ are supported in $\bar{\Omega}_{-}$. The terms with a superscript $s$ are microlocally S waves; and those with a superscript $p$ are P waves.

Denote the restriction of $c_{p}$ and $c_{s}$ to $\bar{\Omega}_{+}$ and $\bar{\Omega}_{-}$, respectively by $c_{p,+}$ or $c_{s,+}$; and $c_{p,-}$, $c_{s,-}$, respectively. A subscript $b$ denotes a boundary value. We know that $\operatorname{WF}(u_{I,b})$ is in the hyperbolic or the mixed region on $T^{*}S$ w.r.t. the speeds $c_{p,+}$ and $c_{s,+}$ assuming non-trivial incoming S and P waves. This may not be true on the negative side, i.e., with respect to the speeds $c_{p,-}$ and $c_{s,-}$ but as we said above, in this section, we are assuming that $\operatorname{WF}(u_{I,b})$ is in the hyperbolic region with respect to them as well.

The transmission conditions $[u]_{\Gamma}=0$, $[Nu]_{\Gamma}=0$ in (1.1) are equivalent to

[TABLE]

where the $\pm$ superscripts indicate that the corresponding operators act in $\Omega_{\pm}$. We will show next that this system is elliptic for recovery of $w_{R,b}$ and $w_{T,b}$ given $w_{I,b}$. In fact, ellipticity is a consequence of the energy preservation. Take the dot product of the two equations above (recall that we work at a fixed point where the metric is transformed to an Euclidean one). We get

[TABLE]

because it can be shown that $(U^{+}_{\text{in}})^{*}M^{+}_{\text{out}}+(M^{+}_{\text{in}})^{*}U^{+}_{\text{out}}=0$ up to smoothing terms. The latter can be proven in the following way. The quadratic form $\langle U_{\textrm{in}}^{+}w_{I,b},M_{\textrm{in}}^{+}w_{I,b}\rangle$ is proportional to the energy flux of $u_{I}$ through $\mathbf{R}\times\Gamma$ as can be shown by integration by parts: we get $2\Re\int_{\mathbf{R}\times\Gamma}\langle u_{t},\Lambda u\rangle$, see, e.g., [29]. Similarly, the other two forms are proportional to energy fluxes, and the signs, after a multiplication by the same constant, are $+,-,+$. Then if $w_{I,b}=0$ (i.e., if (9.2) is homogeneous), the signs of the forms imply the zero solution only. The cancellation equality above reflects the fact that the incoming and the outgoing wave are microlocally separated. We are not going to prove it this way because below we will get a direct confirmation for the principal symbols, which is what we need.

In matrix form, that system is given by

[TABLE]

We compute the principal symbol of the matrix operator applied to $(w_{R,b},w_{T,b})$. As in the previous section, we work at a fixed point where the boundary metric is chosen to be Euclidean. By the invariance under rotations in the $x^{1}x^{2}$ plane, we can perform the computations when $\xi_{2}=0$, as in the previous section. For the principal amplitude of $w$ on $\Gamma$, we will adopt the following notation: $(SH,SV,P)^{T}$, i.e., $P=a^{p}$ in the notation of the previous section, and $a^{s}=(SH,SV)$ is the decomposition of the principal amplitude of the potential (on the boundary) of the S wave $u^{s}=D\times w^{s}$ into shear-horizontal and shear-vertical terms. We use the subscripts $I,R,T$ for the same purpose as above.

The system (9.4) then decouples into a $4\times 4$ one and a $2\times 2$ one. The $4\times 4$ system has the form

[TABLE]

We use the notations $A_{\text{in}}$ and $A+\text{out}$, see (7.7) with plus or minus superscripts depending on which side of $\Gamma$ they are related to. By Lemma 7.1, $(A_{\text{in}},A_{\text{out}})$ is elliptic.

The second system, describing the reflection and the transmission of SH waves, is

[TABLE]

It has a negative determinant, therefore it is elliptic. This decoupling shows that the SH waves do not convert to other modes and reflect and transmit similarly to acoustic waves. We can write (9.6) as

[TABLE]

Multiply those equations to get

[TABLE]

when all $w$’s are real. If they are complex, we can justify this by the equality $\Re(z-w)(\bar{z}+\bar{w})=|z|^{2}-|w|^{2}$. Without going into details, we mention that this is actually an energy equality of the kind (9.3) with $c_{s,\pm}^{2}$ normalization factors since the column vectors of $U$ in (4.5) are not normalized according to the corresponding speed, $\rho_{\pm}$ are volume element factors, the $(\xi_{3,\pm}^{s})^{2}$ factors come from the contribution of an S wave with principal term proportional to $(\xi_{1},0,-\xi_{3}^{s})\times(1,0,0)=\xi_{3}^{s}(0,1,0)$ to $Nu$; and the extra $\xi_{3,\pm}^{s}$ factor accounts for the angle of incidence of reflection/transmission. Equations (9.6) imply that when $SH_{I}\not=0$, we have $SH_{T}\not=0$; and $SH_{R}=0$ when $\mu_{+}\xi_{3,+}^{s}-\mu_{-}\xi_{3,-}^{s}$ which can happen for a fixed $|\xi^{\prime}|$.

We are going back to the system (9.5). We will transform it into a form used in the geophysics literature. Let $\theta_{+}^{p}$, $\theta_{+}^{s}$, $\theta_{-}^{p}$ and $\theta_{-}^{s}$ be the angles between the normal and $\xi^{p}_{R}$, $\xi^{s}_{R}$, $\xi_{T}^{p}$ and $\xi_{T}^{s}$, respectively, see Figure 6. Note that those angles are in $[0,\pi/2)$ and we exclude the zero ones below just to be able to put the equations into the desired form and to compare them with classical results. The singularity at [math] can be resolved by multiplying the corresponding equations by the appropriate sine functions. Then

[TABLE]

and similarly for the other angles.

Divide the first two equations in (9.5) by $\xi_{1}$ and the last two by $\xi_{1}^{2}$, for $\xi_{1}\not=0$, to put the system in the form $A^{\prime}\mathbf{a}=B^{\prime}\mathbf{b}$ with

[TABLE]

and similarly, $B^{\prime}$ is the most right $4\times 2$ block of $A^{\prime}$ with all minus subscripts replaced by plus ones. Here, $\mathbf{a}=(P_{R},SV_{R},P_{T},SV_{T})^{T}$, $\mathbf{b}=(P_{I},SV_{I})^{T}$. The resulting system is the Knott’s equations [16] derived by Knott in 1899 for a flat boundary and constant coefficients. The form here corresponds to [25]. We write them as

[TABLE]

Following [16], we multiply the corresponding sides of the first and the third equations; then do the same thing with the second and the fourth one and add the results to get

[TABLE]

therefore,

[TABLE]

We used here that $\rho_{+}\sin^{2}\theta_{+}^{s}=(\xi_{1}^{2}/\tau^{2})\mu_{+}$ and similarly for the other terms.

As noted by Knott [16], this is an energy equality, stating that the sums of the energy fluxes of the four generated waves, on a principal symbol level, equals that of the incident one. It is also a version of (9.7).

Equation (9.10) implies that the homogeneous system $A^{\prime}\mathbf{a}=0$ has the zero solution only. Therefore, $A^{\prime}$ is elliptic. Explicit formulas for the solution of this system can be found in [1] for the flat constant coefficient case, and those formulas generalize to our case once we make them invariant.

9.2. The general case with incoming waves from both sides

We assume waves coming from both sides, see Figure 7 some of them possibly evanescent, with Dirichlet (and therefore Cauchy) data of their traces on $\Gamma$ in a small neighborhood of some covector in $T^{*}\Gamma$.

We classify the cases by hyperbolic-hyperbolic (HH), hyperbolic-mixed (HM), mixed-mixed (MM), mixed-elliptic (ME) and elliptic-elliptic (EE) according to the location of the wave front of the Cauchy data on the positive/negative side of $\Gamma$.

9.2.1. The hyperbolic-hyperbolic (HH) case

Assume a wave front set in the hyperbolic region on both sides. This is automatically true if on each side, we have both S and P waves. The construction in section 9.1 then generalizes directly. We are going to denote the incoming and the outgoing solutions $w$ on each side by $w_{\text{in}}^{+}$, $w_{\text{out}}^{+}$, $w_{\text{in}}^{-}$, $w_{\text{in}}^{-}$. The transmission conditions (1.1) then take the form

[TABLE]

compare with (9.2) and (7.3). We use the notation in section 7.1 but we put superscripts $+$ and $-$ depending on the side of $\Gamma$ we work on. We use the notation $(P,SV,SH)$ as above for the principal amplitude of $w$ on $\Gamma$, with the corresponding subscripts and the superscripts. Then (9.11) decouples into the following two equations

[TABLE]

and

[TABLE]

compare to (9.5) and(9.6).

The properties of the SH components are similar to those of acoustic waves at an interfaces, see (3.29) there, and the discussion following it. In particular, there is no mode conversion (on principal symbol level at least, which we study).

As above, we can derive the following energy equality:

[TABLE]

For future reference in the case of evanescent modes, we write (9.14) as

[TABLE]

see (9.8). Written this way, (9.15) holds even if the quantities above are not necessarily real; and the proof requires to multiply the first row of (9.12) by the conjugate of the third one and the same for the second and the fourth ones. This is an energy identity, see the paragraph following (9.3). It says that the combined energy flux of all incoming waves on $\Gamma$ (on principal level) equals that of the outgoing ones.

Lemma 9.1**.**

The matrices $(A_{\textrm{\rm in}}^{+},A_{\textrm{\rm out}}^{+})$, $(A_{\textrm{\rm in}}^{-},A_{\textrm{\rm out}}^{-})$, $(A_{\textrm{\rm in}}^{+},A_{\textrm{\rm in}}^{-})$, $(A_{\textrm{\rm out}}^{+},A_{\textrm{\rm out}}^{-})$ are elliptic. Also, system (9.13) is elliptic for $(SH_{\textrm{\rm in}}^{+},SH_{\textrm{\rm out}}^{+})$, and also for $(SH_{\textrm{\rm in}}^{+},SH_{\textrm{\rm in}}^{-})$.

Proof.

The ellipticity of the first two follows from Lemma 7.1. The ellipticity of the next two follows from the energy equality (9.14). The second statement follows from the fact that the corresponding determinants are negative, and positive, respectively. ∎

Note that the ellipticity of $(A_{\textrm{\rm in}}^{+},A_{\textrm{\rm out}}^{+})$ and $(A_{\textrm{\rm in}}^{-},A_{\textrm{\rm out}}^{-})$ holds in the mixed and in the elliptic case as well by the proof of Lemma 7.1.

This has the following implications (without the claim that none of the amplitudes vanishes so far); compare with the discussion following (3.29). Recall that we assume that the Cauchy data on the boundary is in the hyperbolic region with respect to all four speeds.

• (ii)

For every choice of the four incoming waves, there is a unique solution (ellipticity) for the four outgoing ones. Indeed, (9.15) implies unique solution of the homogeneous problem.

• (ii)

An incoming P wave (without any other incoming waves on either side) creates reflected P and S waves and transmitted P and S waves.

• (iii)

The same is true for an incoming S wave.

• (iv)

[Control] For every choice of a principal amplitude of an outgoing transmitted P wave, one can choose incoming S and P waves which would give that pre-assigned transmitted P wave and no (on the principal level) transmitted S wave. The same is true for incoming P waves.

9.2.2. The hyperbolic-mixed (HM) case

Assume the wave front set of the Cauchy data is in the mixed region in $\Omega_{-}$ but still in the hyperbolic one in $\Omega_{+}$. Since we work in the elliptic region for $c_{p,-}$, we will call the principal amplitude of the corresponding microlocal solution $P_{-}$ (no in/out), see Figure 8.

The approach we follow is the same as above — we want to analyze the system (9.11), and for example solve it for all outgoing waves given the incoming ones by proving ellipticity. What changes is that $\xi_{3,-}^{p}$ becomes pure imaginary, see (5.8). One should also change the sign of $\xi_{3,-}^{p}$ in $A_{\text{out}}$ since there are not plus/minus square roots but those entries will be multiplied by zero below. Then (9.12) reduces to

[TABLE]

see (7.11). The energy equality (9.15) reduces to

[TABLE]

see also (3.27). We get that for any choice of the three incoming waves, the resulting system for the three outgoing ones plus $P^{-}$ is elliptic. Indeed, it is enough to show this for the homogeneous system. If all incoming waves vanish, then (9.17) implies $P_{\text{out}}^{+}=SV_{\text{out}}^{+}=SV_{\text{out}}^{-}=0$. Then the only possible non-zero vector in (9.16) is $P^{-}$ but then we can see directly that (9.16) implies $P^{-}=0$. System (9.13) about the SH waves is unaffected by the ellipticity of the $P$ wave. Therefore, constructing the outgoing solution is a well-posed (elliptic) problem.

As far as control from each side is concerned, on the negative one, where $P^{-}$ lies, the Cauchy data is structured; then so is on the positive side. Therefore, the configuration on the positive side cannot be controlled from the negative one. On the other hand, we can create any hyperbolic configuration on the negative side with appropriate waves on the positive one. In particular, if we want $P^{-}=0$, $SV^{-}_{\text{in}}=0$ and $SV^{-}_{\text{out}}\not=0$, we can take the Cauchy data of it and solve (9.16) for the plus amplitudes since on the positive side, we are in the hyperbolic region and the Cauchy problem is elliptic.

Control for SH waves on principal level is the same as in the acoustic case since those waves do not create reflected/transmitted P or SV waves. Since we defined SH/SV waves on principal level only, and the system for the amplitudes is decoupled only on a principal level a priori, the control question needs a further clarification when evanescent P modes are possible. Let us say that we want to create S waves on the negative side with given principal amplitudes $SV_{\text{in}}^{-}$, $SV_{\text{out}}^{-}$, $SH_{\text{in}}^{-}$, $SH_{\text{out}}^{-}$, $P^{-}$. The argument above says that we can chose the principal amplitudes of the waves on the top to make this happen on a principal symbol level, see Figure 8. Then we fix the six waves on the positive side which have those amplitudes as their full ones in those coordinates. For each one, we need to solve, up to infinite order, a transmission, not a control problem, which is well posed. This would possibly create lower order waves on the negative side but it will not change the principal parts. In particular, if we want $SV_{\text{in}}^{-}=SV_{\text{out}}^{-}=SH_{\text{in}}^{-}$ but $SH_{\text{out}}^{-}\not=0$, this step could create lower order $SV_{\text{in}}^{-}$, $SV_{\text{out}}^{-}$, $SH_{\text{in}}^{-}$ waves. This is not a problem since we will need the principal parts later only. We apply the same argument in the cases below.

9.2.3. The mixed-mixed (MM) case

Then on both sides, the S waves are hyperbolic, and $P^{-}$ and $P^{+}$ are evanescent, see Figure 9. In this case, $\xi_{3,\pm}^{p}$ are pure imaginary, see (5.8). Then there is only one evanescent P wave in $\Omega_{-}$ and one in $\Omega_{+}$ and we omit the subscripts “in/out” for them.

As above, we show below that on a principal symbol level, the energy is carried by the S waves only. We also check directly that the homogeneous problem (no incoming waves) has the trivial solution only, including trivial evanescent modes $P^{-}$ and $P^{+}$. Therefore, we still get a well-posed problem for the outgoing solution.

In (9.16), we can formally set $P_{\text{in}}^{+}=0$, $P_{\text{out}}^{+}=P^{+}$ and in the energy equality (9.17), we remove the $P$ amplitudes to get

[TABLE]

and

[TABLE]

with $\xi_{3,\pm}^{p}$ pure imaginary as in (5.8). We will show that (9.18) is elliptic for $SV_{\text{out}}^{-}$, $SV_{\text{out}}^{+}$, $P^{-}$, $P^{-}$, given $SV_{\text{in}}^{-}$, $SV_{\text{in}}^{+}$. As before, it is enough to show that the homogeneous system is uniquely solvable. This follows from Lemma 7.1 or Lemma 9.1 which remain true in the elliptic and the mixed regions.

The SH waves behave as in the acoustic case, see (3.29) and (3.30) and as in the (HM) case.

Control is possible for the SH waves. Let us say that we want to create SH waves on the negative side with prescribed principal amplitudes $SH_{\text{in}}^{-}$, $SH_{\text{out}}^{-}$ and no other waves there. On principal level, we choose $SH_{\text{in}}^{+}$, $SH_{\text{out}}^{+}$ to achieve that. Then, as above, we chose such S waves on the positive side with those principal amplitudes. Solving the direct transmission problem with (hyperbolic only) sources on the positive side, we may get additional waves on the negative ones as shown on Figure 9, left, but they are lower order.

One can also show that control for SV waves on either side is possible from the other one, which would create evanescent $P^{+}$ and $P^{-}$ modes as well. Indeed, to show that given any $SV^{-}_{\text{in}}$, $SV^{-}_{\text{out}}$, we can choose $SV^{+}_{\text{in}}$, $SV^{+}_{\text{out}}$ creating those waves plus the “byproducts” $P^{-}$, $P^{+}$, we need to show that (9.18) is elliptic for $SV^{+}_{\text{in}}$, $SV^{+}_{\text{out}}$, $P^{-}$, $P^{+}$. A direct but tedious computation shows that the determinant of this system equals

[TABLE]

The algebraic structure of this expression implies that this determinant is not identically zero for all $\xi_{1}$ unless none of the coefficients jump at the interface, and we assumed that this could not happen. Therefore, it could be zero for a discrete set of $\xi_{1}$’s only and then we have control.

9.2.4. The hyperbolic-elliptic (HE) case

Assume that both the P and the S waves on the negative side are evanescent but they are hyperbolic on the plus side, see Figure 10. Then we have full reflection on the positive side with respect to all waves. System (9.16) reduces to

[TABLE]

where $P^{-}$ and $SV^{-}$ are evanescent and $\xi_{3,-}^{p}$ and $\xi_{3,-}^{s}$ are pure imaginary as in (5.8) and (5.9). The energy equality takes the form

[TABLE]

This is similar to the hyperbolic case in the Cauchy boundary value problem, see (7). System (9.20) is elliptic for solving for $SV_{\text{out}}^{-}$, $SV_{\text{out}}^{+}$, $P^{-}$, $SV^{-}$ by (9.21) and Lemma 9.1. The SH waves are treated similarly. They experience a full reflection as in the acoustic case.

9.2.5. The mixed-elliptic (ME) case

Assume that only the $SV^{+}$ waves are hyperbolic. Then we have full reflection of the S wave on the positive side with transmitted evanescent $P^{-}$ and $S^{-}$ waves and mode converted $P^{+}$ one on the positive side, see Figure 11. System (9.16) reduces to

[TABLE]

where $P^{-}$, $SV^{-}$ and $P^{+}$ are evanescent and $\xi_{3,\pm}^{p}$ and $\xi_{3,-}^{s}$ are pure imaginary. The energy equality takes the form

[TABLE]

System (9.22) is elliptic for solving for $SV_{\text{out}}^{+}$, $P^{-}$, $SV^{-}$ by (9.23) and Lemma 9.1. The SH waves are treated similarly. They experience a full reflection as in the acoustic case.

9.2.6. The elliptic-elliptic (EE) case. Stoneley waves

We assume now that all waves on both sides are evanescent. Such solutions cannot be created by S or P waves hitting $\Gamma$ but they could be created by boundary sources. We will sketch the construction of such solution known as Stoneley waves first described by R. Stoneley [34] in 1924 in case of flat boundary and constant coefficients, see also [40] for a curved boundary and constant coefficients.

We call the evanescent amplitudes $P^{-}$, $P^{+}$, $SV^{-}$, $SV^{+}$. Then

[TABLE]

Since $\xi_{3,-}^{p}$, $\xi_{3,-}^{s}$, $\xi_{3,+}^{p}$ and $\xi_{3,+}^{s}$ are all pure imaginary, with a positive imaginary part, the matrices above do not really have outgoing properties and the subscript “out” could be omitted. In this region,

[TABLE]

see (7.7). Then $F:=\det(A_{\text{out}}^{+},-A_{\text{out}}^{-})$ is a positively homogeneous function of $(\tau,\xi_{1})$ of order $6$. Writing $F$ as $\xi_{1}^{6}$ times a function $F_{0}$ of $s:=|\tau|/\xi_{1}$ (and the base point $x^{\prime}$), we get that $(A_{\text{out}}^{+},-A_{\text{out}}^{-})$ is elliptic (again, after adjusting the order of the last two rows from 2 to 1) where $F_{0}(s,x^{\prime})\not=0$. Passing to an invariant formulation as in Section 7.2, we can replace $\xi_{1}$ by $|\xi^{\prime}|$; then $s=|\tau|/|\xi^{\prime}|$ with the norm of $\xi^{\prime}$ being the covector one w.r.t. the metric $g$, which in the isotropic case is the boundary metric induced by the Euclidean one. Then $F$ is a homogeneous symbol. Assume that $F_{0}$ has a simple zero for some $s=c_{\text{St}}$ corresponding to the elliptic-elliptic region, i.e., in $s<\min(c_{s,-},c_{s,+})$. Then $F=\left(\tau^{2}-c_{\text{St}}^{2}|\xi^{\prime}|^{2}\right)\!\tilde{F}$ with $\tilde{F}$ elliptic near $\Sigma_{\text{St}}:=\{\tau^{2}=c_{\text{St}}^{2}|\xi^{\prime}|^{2}\}$. Then $(A_{\text{out}}^{+},-A_{\text{out}}^{-})$ is a $\Psi$DO of real principal type (again, the order can be adjusted to be one for all rows) in the sense of [9]. Singularities on $T^{*}\Gamma$ propagate along the null bicharacteristics of the Hamiltonian $H:=\left(\tau^{2}-c_{\text{St}}^{2}|\xi^{\prime}|^{2}\right)$. This is a wave type of Hamiltonian with a wave speed $c_{\text{St}}$ which is slower that the S and the P speeds on either part of $\Gamma$. A well posed problem would be, for example, one with Cauchy data on $\{t=0\}\times\Gamma$.

For every microlocal solution on $\mathbf{R}_{t}\times\Gamma$, we can use its Dirichlet data to extend it to a microlocal solution on both sides of $\Gamma$ as in Section 5.3, see also Section 8.2. Rayleigh waves, can be considered as a limit case of Stoneley waves.

The function $F_{0}$ does have (simple) zero in some cases, at least. Some examples can be found in Stoneley’s original paper [34].

9.3. Summary

We summarize some of the results above as follows.

• (HH)

the hyperbolic-hyperbolic case: we have both P and S waves on either side; each incoming wave creates two reflected and two transmitted (refracted) ones, with mode conversion.

• (HM)

The hyperbolic-mixed case: on one side there are both P and S waves, on the other one, only S waves exists (as solutions propagating singularities); the P wave is evanescent. On other hand, there is total internal reflection of P waves but they can still create transmitted S waves.

• (HE)

The hyperbolic-elliptic case: the S and the P waves on one side are hyperbolic; the S and the P waves on other side are evanescent. Then there is a full reflection from the first side, and the transmitted waves are only evanescent.

• (MM)

The mixed-mixed case: the S waves on both sides are hyperbolic (propagate singularities); the P waves on both sides are evanescent. In particular, an incoming S wave reflects and refracts; and it creates two evanescent P waves on either side by mode conversion.

• (ME)

The mixed-elliptic case: Only the S wave on one side is hyperbolic. In particular, an incoming S wave reflects; and it creates two evanescent P waves on either side by mode conversion and one “reflected” P evanescent one.

• (EE)

The elliptic-elliptic case: All waves are evanescent. Such waves cannot be created by a P or an S wave hitting $\Gamma$ but it could be created by a boundary source. The transmission problem may lose ellipticity and allow for solutions (Stoneley waves) concentrated near $\Gamma$.

9.4. Justification of the parametrix

In the construction above, we work with microlocal solutions which may have singularities but they, and their first derivatives have traces in timelike surfaces. We assume that solutions have wave front set disjoint from bicharacteristics tangential to some of interfaces which can be achieved by choosing the wave front of their Cauchy data disjoint from projections of such directions in $T^{*}\Gamma$. The later set has a zero measure on $S^{*}\partial\Omega$ for $t$ restricted to any fixed finite interval. The construction actually provides an FIO, mapping $f$ to the microlocal outgoing solution $u$ with that boundary data.

To justify the parametrix, we need to subtract it from the actual solution and show that the difference is smooth up to each interface $\Gamma_{i}$. Such a difference $w$ would solve a non-homogeneous problem

[TABLE]

A slightly weaker version of this claim can be proven, which is sufficient for our purposes. We claim that $w$ is $C^{\infty}$ away from $\mathbf{R}\times\Gamma_{j}$ and $\mathbf{R}\times\partial\Omega$, and indeed is conormal at these two in the precise sense that $w\in H_{\mathrm{b},\mathrm{loc}}^{1,\infty}$, meaning $w$ and its first derivatives are in $L^{2}$ locally, and the same remains true if vector fields tangent to $\mathbf{R}\times\Gamma_{j}$ and $\mathbf{R}\times\partial\Omega$ are applied to these iteratively. While this is standard in the scalar case, a proof for (principally) scalar wave equations, for transmission problems, based on quadratic form considerations, showing regularity relative to the quadratic form domain, is given in [6]. This proof uses b-pseudodifferential operators, introduced by Melrose [17], see also [18], and [6] for a brief summary. The simple observation made in [6, Section 4] is that when one has an internal hypersurface, such as $\mathbf{R}\times\Gamma_{j}$, one can treat it as a boundary for this b-analysis by using b-pseudodifferential operators on each half-space (which are manifolds with boundary) with matching normal operators at the common boundary; this was used in [6, Section 4] to prove propagation of singularities in the principally scalar setting. The elastic problem is not principally scalar, which indeed makes the proof of propagation of singularities significantly more difficult using these tools. However, the propagation of global regularity, in the sense that regularity, as measured by $H_{\mathrm{b},\mathrm{loc}}^{1,m}$ (i.e. the space with $m$ b-, or tangential, derivatives relative to $H^{1}_{\mathrm{loc}}$), propagates from $t<0$ to $t\geq 0$ when the right hand side has regularity in $H_{\mathrm{b},\mathrm{loc}}^{-1,m+1}$ (i.e. the space with $m+1$ b-, or tangential, derivatives relative to $H^{-1}_{\mathrm{loc}}$) is straightforward as it does not require microlocalization; slightly modified energy estimates work. This has been carried out in detail by Katsnelson for the elastic wave equation on manifolds with edges in [14, Chapter 11]. The latter are actually more complicated than our setting as the domain of the operator is more delicate, and an essentially identical method of proof works in our case. We also refer to [15] for a brief summary.

We refer to [41] as well, where boundary regularity in the case of constant parameters has been studied.

10. The inverse problem

Assume that there exist two smooth non-positive functions $\mathsf{x}_{s}$ and $\mathsf{x}_{p}$ in $\Omega$ with $\mathrm{d}\mathsf{x}\not=0$, $\mathsf{x}^{-1}(0)=\partial\Omega$, and $\mathsf{x}^{-1}(-j)=\Gamma_{j}$, $j=1,\dots,k$ where $\mathsf{x}$ is either $\mathsf{x}_{s}$ or $\mathsf{x}_{p}$. Assume that the level sets $\mathsf{x}_{s}^{-1}(c)$, $\mathsf{x}_{p}^{-1}(c)$ are strictly convex w.r.t. the speed $c_{s}$, $c_{p}$, respectively, when viewed from $\Gamma_{0}=\partial\Omega$. Of course, we may have just one such function, i.e., $\mathsf{x}_{s}=\mathsf{x}_{p}$ is possible.

Recall that the foliation condition implies non-trapping as noted in [30], for example. In our case, this means that rays in $\Omega_{j}$ not hitting $\Gamma_{j}$ would hit $\Gamma_{j-1}$ both in the future and in the past.

10.1. Recovery in the first layer $\Omega_{1}$

We show first that we can recover $c_{p}$ and $c_{s}$, and then $\rho$, in the first layer $\Omega_{1}$, i.e., between $\partial\Omega$ and $\Gamma_{1}$. In other words, if $\tilde{\rho}$, $\tilde{\mu}$, $\tilde{\nu}$ is another triple of coefficients which have the same piecewise smooth structure with jumps across some $\tilde{\Gamma}_{j}$, producing the same DN map, then they coincide with the non-tilded ones. In the lemmas below, we need solutions with a single incoming singularity (more precisely, with a single radial ray due to the conic nature of the wave front sets) which we can trace until its branches hit $\partial\Omega$ again. We can do this in two ways: first, we can have $f$ in (2.7) with such a single singularity but when we need a specific polarization, we can achieve that by choosing the potential $w$ appropriately, with that singularity. Since the operators $U_{\text{in}}$ and $U_{\text{out}}$, see (5.6) are elliptic in all regions, then the boundary trace of the potentials would have the same wave front sets as the boundary trace of the solution $u$. Or, one can have $\operatorname{WF}(f)$ in a small set by choosing $\operatorname{WF}(w)$ on the boundary small enough and then pass to a limit when $\operatorname{WF}(f)$ shrinks to a single point. Since the arguments based on SH/SV waves require us to trace the leading singularities, i.e., we want to have a well defined order, working with singularities in a small conic set, for example conormal ones, is more convenient. We assume in this section that $g$ is Euclidean since we will need the results of Rachele [21, 23], and Bhattacharyya [2], see Remark 5.1.

Lemma 10.1**.**

Under the convex foliation assumption, $\Lambda$, known for $T\gg 1$ determines uniquely $\Gamma_{1}$, $c_{s}$ and $c_{p}$ in $\Omega_{1}$. If, in addition, $c_{p}\not=2c_{s}$ pointwise in $\Omega_{1}$, then $\rho$ is uniquely determined in $\Omega_{1}$ as well.

Proof.

In this and in the following proof, we consider another triple $\tilde{\rho}$, $\tilde{\mu}$, $\tilde{\nu}$ with the same $\Lambda$, and show that the corresponding quantities, in this case $\Gamma_{1}$ and the three coefficients, coincide. Sometimes, we say that a certain quantity, for example $c_{s}$, is known or can be recovered in some region to indicate that $c_{s}=\tilde{c}_{s}$ there.

First, by [21], we can recover the full jets of $\rho$, $c_{p}$ and $c_{s}$ on $\partial\Omega$. We will recover the speeds $c_{s}$ and $c_{p}$ first. This follows from [30], in any subdomain separated from $\Gamma_{1}$, i.e., for $-1+\varepsilon\leq\mathsf{x}\leq 0$, $\forall\varepsilon\in(0,1)$, with $\mathsf{x}=\mathsf{x}_{s}$ or $\mathsf{x}=\mathsf{x}_{p}$, and it is also Hölder stable there. Indeed, for every unit P or S geodesic connecting boundary points and not intersecting $\Gamma_{1}$, we can construct a microlocal P or S solution in a small neighborhood of that geodesic, extended a bit outside $\Omega$ where the coefficients are extended smoothly as well. Let $f$ be the Dirichlet data of that solution on $\mathbf{R}_{+}\times\partial\Omega$. Then the outgoing solution $\tilde{u}$ having the same Dirichlet data has the same Neumann data as well. Also, the solution will be a P or an S wave, respectively as well, since this property is determined by the trace of $c_{p}$ and $c_{s}$ on $\partial\Omega$, which we recovered, see the end of section 5.2. Therefore, singularities hitting $\partial\Omega$ from inside, will be the same (a singularity hitting $\partial\Omega$ must create singular Cauchy data by the analysis in Section 3.1). So the scattering relations related to $c_{s}$ and $c_{p}$ are the same as those of $\tilde{c}_{s}$ and $\tilde{c}_{p}$ restricted to those geodesics. Note that this argument requires us to know that the corresponding geodesics for the second system do not hit $\tilde{\Gamma}_{1}$. If they do, we would get reflected waves of both kinds (with a possible exception of specific angles which does not change the argument), and we would not get the same Cauchy data. Another way to exclude such rays is to note that they would create singularities of the lens relation near rays tangent to $\tilde{\Gamma}_{1}$.

This proves that $\Omega_{1}\subset\tilde{\Omega}_{1}$, i.e., $\tilde{\Gamma}_{1}$ is below $\Gamma_{1}$, and that $c_{p}=\tilde{c}_{p}$, $c_{s}=\tilde{c}_{s}$ in $\Omega_{1}$. On the other hand, we can swap $\tilde{\Gamma}_{1}$ and $\Gamma_{1}$ in this argument, therefore $\tilde{\Gamma}_{1}=\Gamma_{1}$. Then $c_{s}$ and $c_{p}$ are uniquely determined there. By [2], one can recover $\rho$ in $\Omega_{1}$ as well under the stated condition, therefore then we can recover $\lambda$, $\mu$, too. ∎

Note that here, and in what follows, we have precise control of $T$ which we do not make explicit. Also, local knowledge of $\Lambda$ up to a smoothing operator yields recovery in an appropriate domain of influence, see also [31] for the case of smooth coefficients.

10.2. Recovery in the second layer $\Omega_{2}$

In the next lemmas, we show that we can recover the two speeds in $\Omega_{2}$ under some conditions. The obstruction to the application of the method (but not necessarily to the uniqueness) is existence of totally reflected P and/or S rays on the interior side of $\Gamma_{1}$ for all times (or for long enough, for the case of data on a finite time interval). Since we need rays converging to tangential ones, the microlocal conditions can be described in terms of the sign of the jumps of the speeds at $\Gamma_{1}$.

In what follows, $c|_{\Gamma^{\pm}}$ denotes the limit of $c(x)$ as $x$ approaches $\Gamma$ from the exterior/interior.

Lemma 10.2**.**

Under the assumption in the first sentence of Lemma 10.1, assume additionally that

[TABLE]

*Then $\Gamma_{2}$ and $c_{s}$ are determined uniquely in (the uniquely determined) $\Omega_{2}$. *

We can interpret (10.1) as strict convexity of $\Gamma_{1}$ w.r.t.  $c_{s}$ with a jump since increasing the speeds with depth guarantees strict convexity of the level surfaces. It guarantees no total full reflection of S ways from $\Omega_{2}$ to $\Omega_{1}$. On the other hand, (10.1) implies

[TABLE]

but the only thing we know about $c_{p}|_{\Gamma_{1}^{+}}$ is that it is greater than $c_{s}|_{\Gamma_{1}^{+}}$. In particular, there could be evanescent S to P or P to S transmission from $\Omega_{2}$ to $\Omega_{1}$; or they all could be hyperbolic.

Proof of Lemma 10.2.

Let $x$, $y$ be on $\Gamma_{1}$ connected by a unit speed S geodesic $\gamma_{0}$ staying between $\Gamma_{1}$ and $\Gamma_{2}$. We take an outgoing microlocal solution $u$ concentrated near $\gamma_{0}$, so that $u$ is singular near $x$ when $t$ is near $t_{1}$; and $t=t_{2}$ corresponds to $y$. We choose $u$ to be an SH wave on $\Gamma_{1}$ near $x\in\Gamma_{1}$, see Figure 12. The SH waves behave as acoustic ones on both sides in the (HH), (MH), (HM) and the (MM) cases on principal symbol level, and all those cases are possible. Recall that our convention is to list the top first; in particular, the (MH) configuration is the (HM) one in Section 9.2.2 with the top and the bottom swapped. To create such a wave, we just need to take an S wave coming from $\partial\Omega$ so that its trace on $\Gamma_{1}$ is SH; this can be done by time reversal. On principal level, there will be no other singularities below $\Gamma_{1}$ until that wave hits $\Gamma_{1}$ again. Then that solution will create singular Cauchy data near $y$ and $t$ near $t_{2}$. It is an S wave but not necessarily an SH one at $y$. At least one of the two waves transmitted back to $\Omega_{1}$ would have non-zero principal amplitude if there are two hyperbolic ones, or if there is an S one only, it would be non-zero by the results of the previous section. Then there will be at least one singularity hitting $\partial\Omega$ (which we allow to leave $\Omega$, as above). On the other hand, there might be other waves hitting $\partial\Omega$ at the same place and time coming from waves at $y$ below $\Gamma_{1}$ which can reflect of refract. Since we allow all those waves to leave $\Omega$ freely, they would have different wave front sets or polarizations, and in particular they cannot cancel or alter the singularity of the Cauchy data generated by the waves coming directly from $y$. The simplest way to see that is to do time reversal from the exterior of $\Omega$ back to $\Omega$.

The speeds $c_{s}$ and $c_{p}$ are the same for both systems in $\Omega_{1}$ by Lemma 10.1. We can assume $t_{1}\gg 1$ so that $u$ is smooth on $\partial\Omega$ for $t<\varepsilon$ for some $\varepsilon>0$. Since the solutions constructed above for both systems have the same Cauchy data on $(0,T)\times\partial\Omega$ and we can choose $T\gg 1$, we conclude that the principal part of $u$ on $\Gamma_{1}$ near $t=t_{1}$ is uniquely reconstructed. Note that this argument does not require recovery of $\rho$ in $\Omega_{1}$ since we only need the principal amplitudes and by [21], they do not depend on $\rho$. There might be other singularities on $\Gamma_{1}$ but we can identify $y$ as the first point a singularity comes back to $\Gamma_{1}$, and we can determine the travel time through $\Omega_{2}$ as well. Taking $y\to x$, we can recover the full jet of $c_{s}$ on $\Gamma_{1}^{+}$ by [30, Lemma 2.1]. Since we now know the S metric on $\Gamma_{1}^{-}$, this is enough to recover the lens relation related to $c_{s}$ on $\Gamma_{1}^{-}$, restricted to rays not hitting $\Gamma_{2}$. By [30], this determines $c_{s}$ in $\Omega_{2}$ uniquely, i.e., $c_{s}=\tilde{c}_{s}$ in $\Omega_{2}\cap\tilde{\Omega}_{2}$.

We remark that the magnitudes of the refracted SH waves into $\Omega_{2}$ at $x$ may vary for each of the two systems since we do not know $\rho_{-}:=\rho|_{\Gamma_{1}^{-}}$; see (9.6) where we can write $\mu_{\pm}=\rho_{\pm}c_{s}^{2}$. Their directions however do not depend on $\rho_{-}$ and each one can vanish only for a specific incidence angle (a priori different for each system).

Finally, $\Gamma_{2}=\tilde{\Gamma}_{2}$ since the presence of the interface $\Gamma_{2}$ would create a singularity of the lens relation of the reflected S wave (plus a possible P wave); which would be detected om $\Gamma_{1}$. ∎

Lemma 10.3**.**

Under the assumptions of Lemma 10.2, assume in addition

[TABLE]

*Then $c_{p}$ is uniquely determined in $\Omega_{2}$. *

Remark 10.1**.**

Conditions (10.2) and (10.3) say that there is no total internal reflection of $P\to P$ and $S\to S$ rays from $\Omega_{2}$ to $\Omega_{1}$. On can still have evanescent transmitted $S\to P$ waves from the interior. More precisely, we have the following two generic cases (excluding $c_{s}|_{\Gamma_{1}^{-}}=c_{p}|_{\Gamma_{1}^{+}}$):

[TABLE]

and

[TABLE]

see also (10.2). Evanescent $S\to P$ transmission from the interior happens when (10.4) holds. This is not a problem for the proof since we recovered $c_{s}$ in $\Omega_{2}$ using $SH$ waves. On the other hand, (10.5) implies that all rays from the interior create transmitted rays, i.e., the wave front on $\Gamma_{1}^{+}$ is in the hyperbolic region.

Proof of Lemma 10.3.

We want to use $P$ rays in $\Omega_{2}$ not hitting $\Gamma_{2}$, in particular having Cauchy data on $\Gamma_{1}^{-}$ in the hyperbolic region $c_{p}|_{\Gamma_{1}^{-}}|\xi^{\prime}|<|\tau|$. By (10.3), that Cauchy data will fall in the hyperbolic region on $\Gamma_{1}^{+}$; in other words, we have the (HH) case. We use the control argument in [4] now. In the (HH) case, near $x$ and $t=t_{1}$, we can create an outgoing P wave in $\Omega_{2}$ with no other S or P waves there; in other words, in Figure 6, only $P^{-}_{\text{out}}\not=0$ among all waves on the bottom. Then we extend the four waves on the top until they leave $\Omega$. At $y$, where that wave hits $\Gamma_{1}$ again near $t=t_{2}$, we can apply the same argument to make sure that there are no reflected rays in $\Omega_{2}$, see Figure 13.

By energy preservation, we cannot have zero principal amplitudes of all four rays above $y$. Then by time reversal from $\partial\Omega$, we would know that there is a singularity on $\Gamma_{1}$ at $y$ and $t=t_{2}$, and we would know its wave front set. Note that we do not require knowledge of $\rho$ in $\Omega_{1}$ and $\Omega_{2}$. In principle, the second (tilded) system may have an S wave starting from $x$ at $t=t_{1}$. By the paragraph following Remark 5.1, we must have a non-trivial P wave near $x$ and $y$ (since we have recovered $c_{s}$ already). The P wave arriving at $y$ at $t=t_{2}$ might a priori be due to an S wave from $x$ in $\Omega_{2}$ which has reflected at $\Gamma_{1}$ and mode converted by this is not possible because this would have created a singularity at a moment in the interval $(t_{1},t_{2})$ but we know that such singularity does not exist for either system. Therefore, this recovers the P travel time from $x$ to $y$. Then we recover $c_{p}$ in $\Omega_{2}$ as in the proof of the previous lemma. ∎

Combining those two lemmas, we get the following.

Theorem 10.1**.**

Assume we have two triples of coefficients $\rho$, $\mu$, $\nu$ and $\tilde{\rho}$, $\tilde{\mu}$, $\tilde{\nu}$; and $\Lambda=\tilde{\Lambda}$ with $T\gg 1$. Assume the foliation condition and (10.1) and (10.3) for each one of them. Then $\Gamma_{1}=\tilde{\Gamma}_{1}$, $\Gamma_{2}=\tilde{\Gamma}_{2}$ and $c_{s}=\tilde{c}_{s}$, $c_{p}=\tilde{c}_{p}$ in $\Omega_{1}\cup\Omega_{2}=\tilde{\Omega}_{1}\cup\tilde{\Omega}_{2}$. Also, if $c_{p}\not=2c_{s}$ in $\Omega$, then $\rho=\tilde{\rho}$ in $\Omega_{1}$.

10.3. Recovery of the speeds in the third, etc., layers

This construction can be extended by induction under appropriate conditions:

Theorem 10.2**.**

Assume we have two triples of coefficients $\rho$, $\mu$, $\nu$ and $\tilde{\rho}$, $\tilde{\mu}$, $\tilde{\nu}$. Let

[TABLE]

Assume the foliation condition in $\Omega_{1}\cup\dots\cup\Omega_{k}$. If $\Lambda=\tilde{\Lambda}$ with $T\gg 1$, then $\Gamma_{j}=\tilde{\Gamma}_{j}$, $j=1,\dots,k$ and $c_{s}=\tilde{c}_{s}$, $c_{p}=\tilde{c}_{p}$ in $\Omega_{1}\cup\dots\cup\Omega_{k}$.

Proof.

We show that one can recover the two speeds in $\Omega_{3}$ and then the theorem follows by induction.

We show that we can recover $c_{s}$ there first following the proof of Lemma 10.2. Fix two points $x$ and $y$ on $\Gamma_{2}$. We keep them close enough so that the S geodesic connecting them does not touch $\Gamma_{3}$ (if there is $\Gamma_{3}$, i.e., if $k\geq 3$). We construct a solution $u$ below $\Gamma_{2}$, between $x$ and $y$, of S type (at principal level, as everywhere in this section), see Figure 14. We chose the solution to be SH at $x$ but this is not essential. At $y$, there might be reflection, transmission and mode conversion to evanescent modes. Then near $x$ and at $y$ (and the corresponding times $t_{1}$ and $t_{2}$), the traces of this solution on $\Gamma_{2}$ is in the (HH), (HM) or the (MM) region by (10.6), with the possible exception of finitely many angles giving rise to tangential rays. On the other hand, on principal level, there are only incoming and reflected S waves at $\Gamma_{2}^{+}$ near $(t_{1},x)$, satisfying the transmission conditions, and we can arrange no incoming waves at $x$ from $\Omega_{3}$ by the control argument for SH waves.

We extend those microlocal solutions to $\Omega_{2}$ and $\Omega_{1}$ first as in the proof of Lemma 10.2. We do this starting from $x$ first. On $\Gamma_{1}$, each of the two S branches (meeting at $x$) are in one of the three regions mentioned above excluding directions of measure zero). In each one of those cases, we can choose four or two waves on the top, i.e., in $\Omega_{1}$ which cancel a reflected wave. At $\Gamma_{1}^{-}$, we decompose all S waves into SH and SV ones. The latter can be treated as acoustic waves and can be controlled from the top. The SV ones can be controlled as well as we showed in sections 9.2.1, 9.2.2 and 9.2.3. In Figure 14, for example, point $a$ corresponds to either an (HH) or an (HM) case; and point $b$ corresponds to an (MM) case; so does point $c$. Then we extend all waves in $\Omega_{1}$ to the exterior of $\Omega$, i.e., we let them leave $\Omega$.

At the point $y$, we do the same for the S and the P wave propagating into $\Omega_{2}$. For the $p$ wave (hitting $\Gamma_{1}$ at $d$ in Figure 14), we are at the (HH) zone at $\Gamma_{1}$, and we apply the control argument we used before.

The so constructed microlocal solution vanishes (in this context, that means that it has no leading order singularities) for $t\ll 0$, and by a shifting $t_{1}$, we may assume that this happens for, say, $t<\varepsilon$ with some $\varepsilon>0$ (we need $\varepsilon$ so that we can do a smooth cutoff between $t=0$ and $t=\varepsilon$ and construct an actual solution with the same singularities). Choose $T\gg 1$ so that all outgoing branches starting from $x$ or $y$ reach $\partial\Omega$ before that time.

We are in the situation of Lemma 10.2 now with $\Gamma_{1}$ playing the role of $\partial\Omega$ there with one difference. We have not recovered $\rho$ in $\Omega_{2}$. We claim however that that near the micolrocal solutions along the rays hitting $x$ and $y$, $\tilde{u}$ (corresponding to the second system) has singularities of the same order as $u$. This follows form the following: the Cauchy data on $\mathbf{R}\times\Gamma_{1}^{+}$ and that on $\mathbf{R}\times\Gamma_{1}^{+}$ are related by the transmission conditions (1.1). It follows by (2.16) that they are related by an elliptic operator depending on $\rho$ (recall that we view the three independent coefficients as being $c_{s}$, $c_{p}$ and $\rho$). That dependence makes the refracted and the reflected amplitudes $\rho$ dependent, but it does not change the property of their principal parts being non-zero (except for specific angles). The same conclusion could have been reached by examining thse qualitative behavior of the solution of the microlocal systems, say (9.11) and (9.13) in the (HH) case, as a function of $\rho_{-}$. Therefore, $\tilde{u}$ has leading order singularities in $\Omega_{2}$ along the same rays as $u$ does. By the proof of Lemma 10.2, $y$ is uniquely recovered as the first time the S wave from $x$ hits $\Gamma_{2}$ again. Then the boundary distance function related to $c_{s}$ in $\Omega_{3}$ is uniquely recovered for $x$ close to $y$, which recovers the jet of $c_{s}$ at $\Gamma_{2}^{-}$. Then we know the $c_{s}$ lens relation as well, along rays not touching $\Gamma_{3}$. As in the proofs above, we can detect where $\Gamma_{3}$ is and also recover $c_{s}$ in $\Omega_{3}$.

The recovery of $c_{p}$ in $\Omega_{3}$ goes along the same lines as the proof of Lemma 10.3, using the arguments above. We create a single P wave below $\Gamma_{2}$ connecting $x$ and $y$ and extend it until it reaches $\partial\Omega$ at both sides. At $\Gamma_{2}$, we are in the (HH) case, each ray, extended upwards, will create four new ones. On the upper surfaces, we can have any of the (HH), (HM) and the (MM) cases as above.

We can recover $\Gamma_{3}$ (if exists) by as in the previous lemmas.

The proof for $k>3$ follows by induction. ∎

Remark 10.2**.**

Recovery of $\rho$ in $\Omega_{j}$, $j\geq 2$ seems delicate. The arguments in [23, 2] require the knowledge of the jet of $\rho$ at the boundary up to order three, which is true on $\partial\Omega$ by [21] but proving this on $\Gamma_{j}^{-}$, $j=1,2,\dots$ seems to be not easy.

10.4. Exploiting mode conversion; the PREM model

In the results above, we needed to ensure no total internal reflection of S or P waves or both, from the interior. The mode conversion was not used to obtain information, it was rather a difficulty we had to overcome. Below we show how one can use mode conversion to recover $c_{p}$ when the P waves are totally reflected but the refracted $S$ wave to the exterior is not.

In the Preliminary Reference Earth Model (PREM) [10], in the Upper and the Lower Mantle, the S and the P speeds increase with depth, “on average”, except on a small interval close to the surface. At the boundary of the Lower Mantle and the Outer Core however, the P speed jumps down with depth, hence it does not satisfy (10.3) on that interface. The S speed jumps down to zero, i.e., the Outer Core is believed to be liquid. This violates (10.1) on that interface (and there are no S waves in the Outer Core anyway). Therefore, the P waves in the Outer Core close to tangent ones to their upper boundary are totally reflected (as P waves only) and the results above do not apply for the recovery of $c_{p}$. In this case we can use mode conversion however because PREM shows that those P waves actually produce transmitted (hyperbolic) S waves into the exterior, i.e., condition (10.7) below holds.

An analysis of a solid-liquid model is certainly possible with the methods we develop but it is beyond the scope of this work (see also [7]). We will sketch arguments based on the dynamical system only assuming no S waves below $\Gamma_{1}$ (formally, $c_{s}=0$ there). Those arguments are not a proof since we assume preservation of the microlocal properties in the limit $c_{s}|_{\Omega_{1}}\to 0$.

Assume

[TABLE]

First, we can determine the two speeds in $\Omega_{1}$ per Lemma 10.1.

To recover $c_{p}$ in $\Omega_{2}$, take a P geodesic in $\Omega_{2}$ connecting $x$ and $y$ on $\Gamma_{2}$, so that it does not hit $\Gamma_{2}$; see Figure 15. We can construct a microlocal solution $u$ near it so that it is obtained by an S wave in $\Omega_{1}$ through mode conversion at $\Gamma_{1}$.

To construct such an incoming solution, we can start with such between $x$ and $y$ and time reverse it. Then we take the S branch in $\Omega_{1}$, which on Figure 15 is represented by the dashed most left incoming ray; and let it propagate. There will be a mode conversion in a a neighborhood of $x$, giving use the desired solution. It will have a non-zero principal level energy except possibly for directions of measure zero. There might be a mode converted reflected P wave at $x$ back to $\Omega_{1}$, not shown on Figure 15. If so, we let it propagate and exit $\Omega_{1}$, similarly to the reflected S wave. There will be a reflected P wave, and a transmitted S wave of non-zero principal energy except possibly for angles of measure zero. There might be a P wave propagating from $y$ into $\Omega_{1}$, not shown on the figure. We let them propagate and exit $\Omega_{1}$ through $\partial\Omega$.

Since the tilded system has the same Cauchy data on $(0,T)\times\partial\Omega$ (as above, we shift the time, if needed so that the solution is smooth for $t<0$), and $c_{p}=\tilde{c}_{p}$, $c_{s}=\tilde{c}_{2}$ in $\Omega_{1}$ by Lemma 10.1, we get that the principal part of $u$ and $\tilde{u}$ coincide in the domain of influence. We recall that the principal parts do not depend on $\rho$. Then $u$ and $\tilde{u}$ have the same Cauchy data on $\Gamma_{1}^{+}$ as well. Then we can identify $y$ by the point where the first (in time) singularity hits $\Gamma_{1}$ again. The rest is as in the proof of the previous results.