Determination of the density in a nonlinear elastic wave equation

Gunther Uhlmann Department of Mathematics, University of Washington, Seattle, WA 98195, USA; Institute for Advanced Study, The Hong Kong University of Science and Technology, Kowloon, Hong Kong, China ([email protected]) and Jian Zhai School of Mathematical Sciences, Fudan University, Shanghai 200433, China ([email protected]).

Abstract.

This is a continuation of our study [41] on an inverse boundary value problem for a nonlinear elastic wave equation. We prove that all the linear and nonlinear coefficients can be recovered from the displacement-to-traction map, including the density, under some natural geometric conditions on the wavespeeds.

Key words and phrases:

elastic waves, inverse boundary value problem, quasilinear equation, Gaussian beams

J. Zhai is supported by National Key Research and Development Programs of China (No. 2023YFA1009103), Science and Technology Commission of Shanghai Municipality (23JC1400501)

1. Introduction

Let $\Omega\subset\mathbb{R}^{3}$ be a bounded domain with smooth boundary $\partial\Omega$ . Denote $x=(x_{1},x_{2},x_{3})$ to be the Cartesian coordinates. Here $\Omega$ represents an elastic, nonhomogeneous, isotropic object with Lamé parameters $\lambda,\mu$ and density $\rho$ . We take into account the nonlinear behavior of elastic medium in this article. Let the vector $u=(u_{1},u_{2},u_{3})$ be the displacement vector, and the (nonlinear) strain tensor is

\varepsilon_{ij}(u)=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}+\frac{\partial u_{k}}{\partial x_{i}}\frac{\partial u_{k}}{\partial x_{j}}\right).

We consider the nonlinear model of elasticity as in [11], whereas the stress tensor is of the form,

(1)

\begin{split}S_{ij}=&\lambda\varepsilon_{mm}\delta_{ij}+\lambda\widetilde{\varepsilon}_{mm}\frac{\partial u_{i}}{\partial x_{j}}+2\mu\left(\varepsilon_{ij}+\widetilde{\varepsilon}_{jn}\frac{\partial u_{i}}{\partial x_{n}}\right)\\ &+\mathscr{A}\widetilde{\varepsilon}_{in}\widetilde{\varepsilon}_{jn}+\mathscr{B}(2\widetilde{\varepsilon}_{nn}\widetilde{\varepsilon}_{ij}+\widetilde{\varepsilon}_{mn}\widetilde{\varepsilon}_{mn}\delta_{ij})+\mathscr{C}\widetilde{\varepsilon}_{mm}\widetilde{\varepsilon}_{nn}\delta_{ij}+\mathcal{O}(u^{3}),\end{split}

where $\widetilde{\varepsilon}$ is the linearized strain tensor

\widetilde{\varepsilon}_{ij}(u)=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right),

For the rest of the paper, we will just discard the $\mathcal{O}(u^{3})$ terms in (1).

We consider the initial boundary value problem

(2)

\begin{split}&\rho\frac{\partial^{2}u}{\partial t^{2}}-\nabla\cdot S(x,u)=0,\quad(t,x)\in(0,T)\times\Omega,\\ &u(t,x)=f(t,x),\quad(t,x)\in(0,T)\times\partial\Omega,\\ &u(0,x)=\frac{\partial}{\partial t}u(0,x)=0,\quad x\in\Omega,\end{split}

where $S$ is of the form (1). Denote the displacement-to-traction map as

\Lambda:f=u|_{(0,T)\times\partial\Omega}\rightarrow\nu\cdot S(x,u)|_{(0,T)\times\partial\Omega},

where $\nu$ is the exterior normal unit vector to $\partial\Omega$ . We consider the inverse problem of determining $\lambda,\mu,\rho,\mathscr{A},\mathscr{B},\mathscr{C}$ from $\Lambda$ . The well-posedness of the Dirichlet problem with small boundary data is established in [11].

Closely related is the inverse boundary value problem for the linear elastic wave equations

(3)

\begin{split}&\rho\frac{\partial^{2}u}{\partial t^{2}}-\nabla\cdot S^{L}(x,u)=0,\quad(t,x)\in(0,T)\times\Omega,\\ &u(t,x)=f(t,x),\quad(t,x)\in(0,T)\times\partial\Omega,\\ &u(0,x)=\frac{\partial}{\partial t}u(0,x)=0,\quad x\in\Omega,\end{split}

for which one wants to recover $\lambda,\mu,\rho$ from the (linearized) Dirichlet-to-Neumann (DtN) map defined as

\Lambda^{lin}:f=u|_{(0,T)\times\partial\Omega}\rightarrow\nu\cdot S^{L}(x,u)|_{(0,T)\times\partial\Omega}.

Here $S^{L}$ is the linearized stress

S^{L}_{ij}(x,u)=\lambda\widetilde{\varepsilon}_{mm}\delta_{ij}+2\mu\widetilde{\varepsilon}_{ij}.

Throughout the paper, we assume the strong ellipticity condition

(4)

\mu>0,\quad 3\lambda+2\mu>0~{}\text{on}~{}\overline{\Omega}.

Denote

c_{P}=\sqrt{\frac{\lambda+2\mu}{\rho}},\quad c_{S}=\sqrt{\frac{\mu}{\rho}},

which are actually the wavespeeds for P- and S- waves respectively. By the assumption (4), we have $c_{P}>c_{S}$ in $\overline{\Omega}$ . Denote the Riemannian metrics associated with the $P/S$ - wavespeeds to be

g_{P/S}=c_{P/S}^{-2}\mathrm{d}s^{2},

where $\mathrm{d}s^{2}$ is the Euclidean metric. Now, one can consider $(\Omega,g_{\bullet})$ as a compact Riemannian manifold with boundary $\partial\Omega$ , where $\bullet=P,S$ .

We assume that $(\Omega,g_{\bullet})$ is non-trapping and $\partial\Omega$ is convex with respect to the metric $g_{\bullet}$ . Denote $\mathrm{diam}_{P/S}(\Omega)$ be the diameter of $\Omega$ with respect to $g_{P/S}$ . More precisely,

\mathrm{diam}_{P/S}(\Omega)=\sup\{\textit{lengths of all geodesics in }(\Omega,g_{P/S})\}.

The main result of this paper is the following theorem.

Theorem 1.

Assume $\partial\Omega$ is strictly convex with respect to $g_{P/S}$ , and either of the following conditions holds

(1)

$(\Omega,g_{P/S})$ is simple;
(2)

$(\Omega,g_{P/S})$ satisfies the foliation condition.

If $T>2\,\max\{\mathrm{diam}_{S}(\Omega),\mathrm{diam}_{P}(\Omega)\}$ , then $\Lambda$ determines $\lambda,\mu,\rho,\mathscr{A},\mathscr{B},\mathscr{C}$ in $\overline{\Omega}$ uniquely.

Assume $(\Omega,g)$ is a compact Riemannian manifold with strictly convex boundary $\partial\Omega$ . We recall that (1) $(\Omega,g)$ called to be simple if any two points $x,y\in\Omega$ can be connected by a unique geodesic, contained in $\Omega$ , depending smoothly on $x$ and $y$ . (2) $(\Omega,g)$ satisfies the foliation condition if there is a smooth strictly convex function $f:M\rightarrow\mathbb{R}$ . For more discussions on the foliation condition, we refer to [27]. We emphasize here that either of the two geometrical conditions implies that $(\Omega,g)$ is non-trapping. The foliation condition is satisfies if the manifold has negative sectional curvatures, or non-positive sectional curvatures if the manifolds is simply connected. The foliation condition allows for conjugate points.

For the inverse boundary value problem for the linear elastic wave equation (3), the uniqueness of the wavespeeds $c_{P/S}$ , under the above two geometrical conditions, is proved in [29, 35] respectively. The problem is reduced to the lens rigidity problem for a Riemannian metric which is conformally Euclidean. The related lens rigidity problem is solved under the simplicity condition in [25] and under the foliation condition in [34]. To recover $\lambda,\mu,\rho$ simultaneously, the uniqueness of $\rho$ needs to be proved additionally. For the uniqueness of the density, the problem can be reduced to the geodesic ray transform of a 2-tensor using P-wave measurements (cf. [30, 6]), but this reduction requires a technical assumption $c_{P}\neq 2c_{S}$ . The $s$ -injectivity of this geodesic ray transform is proved under the strictly convex foliation condition [36]. Under simplicity condition, the $s$ -injectivity is proved under extra curvature conditions [33, 10, 26].

In our previous work [41], we proved the uniqueness of $\mathscr{A},\mathscr{B},\mathscr{C}$ assuming $\lambda,\mu,\rho$ are known. For the uniqueness of $\rho$ , extra conditions are needed if one resorts to the result for the linear equation as mentioned above. In this paper, we will utilize the nonlinearity to show the uniqueness of $\rho$ , without any additional assumptions.

The proof is based on the construction of Gaussian beam solutions. Compared to our previous work [41], since now the density is not known, the reflection of waves at the boundary is hard to control. We remark here that elastic waves would undergo a $P/S$ mode conversion when reflected. We refer to the work [37] for this reflection behavior. The reflection of Gaussian beam solutions is carefully characterized below in Section 3.

The nonlinear interaction of distorted planes was first used in [21] to recover a Lorentzian metric in a semilinear wave equation. Since then, many inverse problems for nonlinear hyperbolic equations have been studied. See [40] for a review. We also refer to [24, 11, 39, 15, 41, 20, 43, 8, 7, 32, 23, 22, 2, 4, 31, 17, 16, 42, 9] and the references therein.

The rest of the paper is organized as follows. In Section 2, we construct Gaussian beam solutions for the linear elastic wave equations. In Section 3, the reflection of Gaussian beams on the boundary is carefully characterized. In Section 4, we give the proof of the main result.

2. Gaussian beams

In this section, we construct Gaussian beam solutions to the linear elastic wave equation

\mathcal{L}_{\lambda,\mu,\rho}u:=\rho\frac{\partial^{2}u}{\partial t^{2}}-\nabla\cdot S^{L}(x,u)=0,

Much of the work has been carried out in [41]. Gaussian beams have also been used to study various inverse problems for both elliptic and hyperbolic equations [19, 3, 5, 12, 13, 15, 14]. For the construction, we need to introduce the Fermi coordinates.

2.1. Fermi coordinates

Denote

M=\mathbb{R}\times\Omega,

and consider $M$ as a Lorentzian manifold with metric $\overline{g}_{\bullet}:=-\mathrm{d}t^{2}+g_{\bullet}=-\mathrm{d}t^{2}+c_{\bullet}^{-2}\mathrm{d}s^{2}$ , where $\bullet=P$ or $S$ . We can extend $g_{\bullet}$ smoothly to a slightly larger domain $\Omega^{\prime}$ , such that $\Omega\subset\subset\Omega^{\prime}$ , and consider the Lorentzian manifold $(M^{\prime},\overline{g}_{\bullet}):=(\mathbb{R}\times\Omega^{\prime},\overline{g}_{\bullet})$ .

Consider a null-geodesic $\vartheta$ in $(M^{\prime},\overline{g}_{\bullet})$ . Notice that $\vartheta$ can be expressed as $\vartheta(t)=(t,\gamma(t))$ , where $\gamma$ is a unit-speed geodesic in the Riemannian manifold $(\Omega^{\prime},g_{\bullet})$ . Assume that $\vartheta$ passes through two boundary points $(t_{-},\gamma(t_{-}))$ and $(t_{+},\gamma(t_{+}))$ of $\partial M$ , where $t_{-},t_{+}\in(0,T)$ , $t_{-}<t_{+}$ and $\gamma(t_{-}),\gamma(t_{+})\in\partial\Omega$ . We assume $\vartheta:[t_{-}-\epsilon,t_{+}+\epsilon]\rightarrow M^{\prime}$ , by extending it if necessary, and introduce the Fermi coordinates in a neighborhood of $\vartheta$ . We will follow the construction of the coordinates in [13]. See also [20], [39].

Assume $\gamma(t_{0})=x_{0}\in\Omega$ where $t_{0}\in(t_{-},t_{+})$ . Choose $\alpha_{2},\alpha_{3}\in T_{x_{0}}\Omega^{\prime}$ such that $\{\dot{\gamma}(t_{0}),\alpha_{2},\alpha_{3}\}$ forms an orthonormal basis for $T_{x_{0}}\Omega^{\prime}$ . Let $s$ denote the arc length along $\gamma$ from $x_{0}$ . We note here that $s$ can be positive or negative, and $(t_{0}+s,\gamma(t_{0}+s))=\vartheta(t_{0}+s)$ . For $k=2,3$ , let $e_{k}(s)\in T_{\gamma(t_{0}+s)}\Omega^{\prime}$ be the parallel transport of the vector $\alpha_{k}$ along $\gamma$ to the point $\gamma(t_{0}+s)$ .

Define the coordinate system $(y^{0}=t,y^{1}=s,y^{2},y^{3})$ through $\mathcal{F}_{1}:\mathbb{R}^{1+3}\rightarrow\mathbb{R}\times\Omega^{\prime}$ :

\mathcal{F}_{1}(y^{0}=t,y^{1}=s,y^{2},y^{3})=(t,\exp_{\gamma(t_{0}+s)}\left(y^{2}e_{2}(s)+y^{3}e_{3}(s)\right)).

Then the null-geodesic $\vartheta$ can be expressed as

\vartheta=\{t-s=t_{0},y^{2}=y^{3}=0\}.

Notice that $(y^{1}=s,y^{2},y^{3})$ gives a coordinate system on $\Omega^{\prime}$ in a neighborhood of $\gamma$ such that

g_{\bullet}|_{\gamma}=\sum_{j=1}^{3}(\mathrm{d}y^{j})^{2},\quad\text{and}\quad\frac{\partial g_{\bullet,jk}}{\partial y^{i}}\Big{|}_{\gamma}=0,~{}1\leq i,j,k\leq 3.

Under this coordinate system, the Euclidean metric $g_{E}:=\mathrm{d}s^{2}$ in a neighborhood of $\gamma$ takes the form

g_{E}=c_{\bullet}^{2}g_{\bullet}=\sum_{1\leq i,j\leq 3}c^{2}_{\bullet}g_{\bullet,ij}\mathrm{d}y^{i}\mathrm{d}y^{j},

and the Christoffel symbols for the Euclidean metric are given by

(5)

\begin{split}&\Gamma_{\alpha\beta}^{1}=-c_{\bullet}^{-1}\frac{\partial c_{\bullet}}{\partial s}g_{\bullet,\alpha\beta},\quad\Gamma_{1\alpha}^{\beta}=\delta^{\alpha}_{\beta}c_{\bullet}^{-1}\frac{\partial c_{\bullet}}{\partial s},\quad\Gamma_{1\alpha}^{1}=c_{\bullet}^{-1}\frac{\partial c_{\bullet}}{\partial y^{\alpha}},\\ &\Gamma_{11}^{\alpha}=-c_{\bullet}^{-1}g_{\bullet}^{\alpha\beta}\frac{\partial c_{\bullet}}{\partial y^{\beta}},\quad\Gamma_{11}^{1}=c_{\bullet}^{-1}\frac{\partial c_{\bullet}}{\partial s},\end{split}

where $\alpha,\beta\in\{2,3\}$ .

Introduce the map $(z^{0},z^{\prime}):=(z^{0},z^{1},z^{2},z^{3})=\mathcal{F}_{2}(y^{0}=t,y^{1}=s,y^{2},y^{3})$ , where

z^{0}=\tau=\frac{1}{\sqrt{2}}(t-t_{0}+s),\quad z^{1}=r=\frac{1}{\sqrt{2}}(-t+t_{0}+s),\quad z^{j}=y^{j},\,j=2,3.

The Fermi coordinates $(z^{0}=\tau,z^{1}=r,z^{2},z^{3})$ near $\vartheta$ is given by $\mathcal{F}:\mathbb{R}^{1+3}\rightarrow\mathbb{R}\times\Omega^{\prime}$ , where $\mathcal{F}=\mathcal{F}_{1}\circ\mathcal{F}_{2}^{-1}$ . Denote $\tau_{\pm}=\sqrt{2}(t_{\pm}-t_{0})$ . Then on $\vartheta$ we have

\overline{g}_{\bullet}|_{\vartheta}=2\mathrm{d}\tau\mathrm{d}r+\sum_{j=2}^{3}(\mathrm{d}z^{j})^{2}\quad\text{and}\quad\frac{\partial\overline{g}_{\bullet,jk}}{\partial z^{i}}\Big{|}_{\vartheta}=0,~{}0\leq i,j,k\leq 3.

2.2. Gaussian beams

With $\varrho$ as a large parameter, the Gaussian beam solutions have the asymptotic form as

u_{\varrho}=\mathfrak{a}e^{\mathrm{i}\varrho\varphi},

with

\varphi=\sum_{k=0}^{N}\varphi_{k}(\tau,z^{\prime}),\quad\mathfrak{a}(\tau,z^{\prime})=\chi\left(\frac{|z^{\prime}|}{\delta}\right)\sum_{k=0}^{N+1}\varrho^{-k}\mathbf{a}_{k}(\tau,z^{\prime}),\quad\mathbf{a}_{k}(\tau,z^{\prime})=\sum_{j=0}^{N}\mathbf{a}_{k,j}(\tau,z^{\prime}),

which is compactly supported in a neighborhood of $\vartheta$ ,

\mathcal{V}=\{(\tau,z^{\prime})\in M^{\prime}|\tau\in\left[\tau_{-}-\frac{\epsilon}{\sqrt{2}},\tau_{+}+\frac{\epsilon}{\sqrt{2}}\right],\,|z^{\prime}|<\delta\}.

Here $\delta>0$ is a small parameter. For each $j$ , $\varphi_{j}$ and $\mathbf{a}_{k,j}$ are complex valued homogeneous polynomials of degree $j$ with respect to the variables $z^{i}$ , $i=1,2,3$ . The smooth function $\chi:\mathbb{R}\rightarrow[0,+\infty)$ satisfies $\chi(t)=1$ for $|t|\leq\frac{1}{4}$ and $\chi(t)=0$ for $|t|\geq\frac{1}{2}$ . We refer to [15] for more details. Under the coordinates $(y^{1},y^{2},y^{3})$ on $\Omega^{\prime}$ , we denote the (co)vector $\mathbf{a}_{k}=(a_{k1},a_{k2},a_{k3})$ . The parameter $\delta$ is small enough such that $\mathfrak{a}|_{t=0}=\mathfrak{a}|_{t=T}=0$ .

In the following calculations, we work under the coordinate system $(y^{1},y^{2},y^{3})$ on the Riemannian manifold $(\Omega^{\prime},g_{E})$ . That is, all the inner products and covariant derivatives are with respect to the Euclidean metric $g_{E}$ . For simplicity of notations, we denote $\overline{g}=\overline{g}_{\bullet}$ , $g=g_{\bullet}$ and $c=c_{\bullet}$ , $\bullet=S,P$ . In a neighborhood of $\vartheta$ , we can write (cf. [41])

(6)

\rho\partial_{t}^{2}u_{\varrho}-\nabla\cdot S^{L}(u_{\varrho})=e^{\mathrm{i}\varrho\varphi}\left(\varrho^{2}\mathcal{I}_{1}+\sum_{k=0}^{N}\varrho^{1-k}\mathcal{I}_{k+2}+\mathcal{O}(\varrho^{-N})\right),

where

\mathcal{I}_{1}=-\rho(\partial_{t}\varphi)^{2}\mathbf{a}_{0}+(\lambda+\mu)\langle\mathbf{a}_{0},\nabla\varphi\rangle\nabla\varphi+\mu|\nabla\varphi|^{2}\mathbf{a}_{0},

and

\begin{split}(\mathcal{I}_{2})_{i}=&\rho(\partial^{2}_{t}\varphi)a_{0i}+2\rho\partial_{t}\varphi\partial_{t}a_{0i}+\mathrm{i}\rho(\partial_{t}\varphi)^{2}a_{1i}\\ &-c^{-2}\left(\partial_{i}(\lambda a_{0k}\varphi_{;\ell}g^{k\ell})+\lambda a_{0k}\varphi_{;\ell}g^{k\ell}\partial_{m}g_{ij}g^{jm}+\partial_{m}(\mu a_{0i}\varphi_{;j}+\mu a_{0j}\varphi_{;i})g^{jm}\right)\\ &-c^{-2}\left(\lambda a_{0k;\ell}g^{k\ell}\varphi_{;i}+\mu(a_{0i;j}+a_{j;i})\varphi_{;m}g^{jm}+\mathrm{i}\lambda a_{1k}\varphi_{;\ell}g^{k\ell}\varphi_{;i}+\mathrm{i}\mu(a_{1i}\varphi_{;j}\varphi_{;m}g^{jm}+\varphi_{;i}a_{1j}\varphi_{;m}g^{jm})\right)\\ &+\Gamma^{n}_{im}g^{mj}c^{-2}(\lambda a_{0k}\varphi_{;\ell}g^{k\ell}g_{nj}+\mu a_{0n}\varphi_{;j}+\mu a_{0j}\varphi_{;n})\\ &+\Gamma^{n}_{jm}g^{mj}c^{-2}(\lambda a_{0k}\varphi_{;\ell}g^{k\ell}g_{ni}+\mu a_{0n}\varphi_{;i}+\mu a_{0i}\varphi_{;n}).\end{split}

We will need to construct the phase function $\varphi$ and the amplitude $\mathfrak{a}$ such that

(7)

\frac{\partial^{\Theta}}{\partial z^{\Theta}}\mathcal{I}_{k}=0\text{ on }\vartheta

for $\Theta=(0,\Theta_{1},\Theta_{2},\Theta_{3})$ with $|\Theta|\leq N$ and $k=1,2,\cdots,N+2$ .

2.3. S-waves

For the construction of S-waves, we note that $g=g_{S}$ . In order for $\mathcal{I}_{1}$ to vanish up to order $N$ on $\vartheta$ (cf. (7) with $k=1$ ), we take $\varphi$ and $\mathbf{a}_{0}$ such that

(8)

\frac{\partial^{\Theta}}{\partial y^{\Theta}}(\mathcal{S}\varphi)=0\text{ on }\vartheta,

where $\mathcal{S}\varphi=\mu|\nabla\varphi|^{2}-\rho(\partial_{t}\varphi)^{2}$ , and

(9)

\frac{\partial^{\Theta}}{\partial y^{\Theta}}\langle\mathbf{a}_{0},\nabla\varphi\rangle=0\text{ on }\vartheta,

for any $|\Theta|\leq N$ . Taking $\Theta=0$ in (9), we conclude that $a_{01}|_{\vartheta}=0$ .

For the construction of the phase function $\varphi$ , we can take

\varphi=\sum_{k=0}^{N}\varphi_{k}(\tau,z^{\prime})

such that

\varphi_{0}(\tau,r,z^{\prime\prime})=0,\quad\varphi_{1}(\tau,r,z^{\prime\prime})=r,\quad\varphi_{2}(\tau,z^{\prime})=\sum_{i,j=1}^{3}H_{ij}(\tau)z^{i}z^{j},

where $H$ is a symmetric matrix solving the Riccati equation

(10)

\frac{\mathrm{d}}{\mathrm{d}\tau}H+HCH+D=0,\tau\in\left(\tau_{-}-\frac{\epsilon}{2},\tau_{+}+\frac{\epsilon}{2}\right),\quad H(\tau_{0})=H_{0},\text{ with }\Im H_{0}>0,

where $C$ , $D$ are matrices with $C_{11}=0$ , $C_{ii}=2$ , $i=2,3$ , $C_{ij}=0$ , $i\neq j$ , $D_{ij}=\frac{1}{4}(\partial_{ij}^{2}\overline{g}^{11})$ , $H_{0}$ is any given symmetric matrix with $\Im H_{0}>0$ . Here $(\overline{g}^{ij})$ is the inverse of the Lorentzian metric $(\overline{g}_{ij})$ under the Fermi coordinates $(z^{0},z^{1},z^{2},z^{3})$ , i.e., $\overline{g}_{ij}=\overline{g}(\frac{\partial}{\partial z^{i}},\frac{\partial}{\partial z^{j}})$ . Then the equation (10) has a unique solution with $\Im(H(\tau))>0$ for all $\tau$ .

For solving the Ricatti equation (10), we take

H(\tau)=Z(\tau)Y(\tau)^{-1},

where $Z(\tau)$ and $Y(\tau)$ are solutions to the first order linear ODEs

\begin{split}\frac{\mathrm{d}}{\mathrm{d}\tau}Y=CZ,\quad Y(\tau_{0})=Y_{0},\\ \frac{\mathrm{d}}{\mathrm{d}\tau}Z=-DY,\quad Z(\tau_{0})=H_{0}Y_{0},\end{split}

where $Y_{0}$ is any non-degenerate matrix. Here $Y(\tau)$ is non-degnerate for all $\tau$ . Moreover, the following identity holds

\det(\Im H(\tau))|\det(Y(\tau))|^{2}=c_{0},

where $c_{0}$ is a constant independent of $\tau$ . For more discussions on the Ricatti equation, we refer to [18]. The higher order terms $\varphi_{k}$ , $k=3,\cdots,N,$ can be constructed so that (8) is satisfied (cf. [13, 15]).

Next consider the equation $(\mathcal{I}_{2})_{i}=0$ for $i=2,3$ . We have, for $\alpha=2,3$ ,

\begin{split}\left(\rho\partial^{2}_{t}\varphi-c_{S}^{-2}\mu\partial_{s}\varphi-c_{S}^{-2}\mu\sum_{\beta=2}^{3}\frac{\partial^{2}\varphi}{\partial y^{\beta}\partial y^{\beta}}\right)a_{0\alpha}\\ -\sqrt{2}\left(\rho\partial_{t}a_{0\alpha}+c_{S}^{-2}\mu\frac{\partial_{s}a_{0\alpha}}{\partial s}\right)&\\ -\frac{1}{\sqrt{2}}c_{S}^{-2}\partial_{s}\mu a_{0\alpha}+\frac{1}{\sqrt{2}}\mu c_{S}^{-3}\frac{\partial c_{S}}{\partial s}a_{0\alpha}+\frac{1}{2}\mathrm{i}a_{1\alpha}(\rho-c_{S}^{-2}\mu)&\\ -c_{S}^{-2}(\lambda+\mu)\left(\frac{1}{\sqrt{2}}\frac{\partial a_{01}}{\partial y^{\alpha}}+\sum_{\beta=2}^{3}a_{0\beta}\frac{\partial^{2}\varphi}{\partial y^{\alpha}\partial y^{\beta}}\right)&=0.\end{split}

Notice that the above equation can be rewritten as

(11)

\mathcal{T}a_{0\alpha}+\frac{1}{2}\mathrm{i}a_{1\alpha}(\rho-c_{S}^{-2}\mu)-c_{S}^{-2}(\lambda+\mu)\left(\frac{1}{\sqrt{2}}\frac{\partial a_{01}}{\partial y^{\alpha}}+\sum_{\beta=2}^{3}a_{0\beta}\frac{\partial^{2}\varphi}{\partial y^{\alpha}\partial y^{\beta}}\right)=0,

where

\mathcal{T}=2\frac{\partial}{\partial\tau}+\left[\frac{1}{\mu}\frac{\partial\mu}{\partial\tau}-c_{S}^{-1}\frac{\partial c_{S}}{\partial\tau}+\frac{1}{\det(Y_{S})}\frac{\partial\det(Y_{S})}{\partial\tau}\right].

Consider the equation (9) with $|\Theta|=1$ , we have

(12)

\frac{1}{\sqrt{2}}\frac{\partial a_{01}}{\partial y^{k}}+\sum_{\beta=2}^{3}a_{0\beta}\frac{\partial^{2}\varphi}{\partial y^{k}\partial y^{\beta}}=0,

for $k=1,2,3$ . Insert (12) (with $k=2,3$ ) into (11), together with the fact $\rho=c_{S}^{-2}\mu$ , we have the following transport equation for $a_{0\alpha}$ , $\alpha=2,3$ ,

\mathcal{T}a_{0\alpha}=0.

By solving the above transport equation, one can determine $a_{0\alpha}|_{\vartheta}$ for $\alpha=2,3$ . In particular we can take

a_{0\alpha}(\tau)=c_{\alpha}\det(Y_{S}(\tau))^{-1/2}c_{S}(\tau,0)^{-1/2}\rho(\tau,0)^{-1/2},

where the constant $c_{\alpha}$ depends on the initial value of $a_{0\alpha}$ . Using the equation (12) again, we can now determine $\frac{\partial^{\Theta}a_{01}}{\partial y^{\Theta}}\Big{|}_{\vartheta}$ for $|\Theta|=1$ .

Noting that we have determined $\mathbf{a}_{0}$ and $\partial a_{01}$ on $\vartheta$ and using the equation $(\mathcal{I}_{2})_{i=1}=0$ , we end up with

\begin{split}\mathrm{i}\rho(\partial_{t}\varphi)^{2}a_{11}-\mathrm{i}c_{S}^{-2}(\lambda a_{1k}\varphi_{;\ell}g^{k\ell}\varphi_{;1}+\mu a_{11}\varphi_{;j}\varphi_{;m}g^{jm}+\varphi_{;1}a_{1j}\varphi_{;m}g^{jm})&\\ +2\rho\partial_{t}\varphi\partial_{t}a_{01}-c_{S}^{-2}(\lambda\partial_{1}a_{0k}\varphi_{;\ell}g^{k\ell}+(\mu\partial_{m}a_{01}\varphi_{;j}+\mu\partial_{m}a_{0j}\varphi_{;1})g^{jm})&\\ -c_{S}^{-2}(\lambda\partial_{\ell}a_{0k}g^{k\ell}\varphi_{;1}+\mu(\partial_{j}a_{01}+\partial_{1}a_{0j})\varphi_{;m}g^{jm})&=F(\mathbf{a}_{0}),\end{split}

which further simplifies to

(13)

\mathrm{i}(\rho-c_{S}^{-2}(\lambda+2\mu))a_{11}-\sqrt{2}c_{S}^{-2}(\lambda+\mu)(\partial_{2}a_{02}+\partial_{3}a_{03})=F(\mathbf{a}_{0},\partial a_{01}).

Notice that $\rho-c_{S}^{-2}(\lambda+2\mu)\neq 0$ . We will also use equation (9) with $|\Theta|=2$ , which can be written as

(14)

\frac{1}{\sqrt{2}}\frac{\partial^{2}a_{01}}{\partial y^{m}\partial y^{k}}+\frac{\partial a_{0j}}{\partial y^{k}}\frac{\partial^{2}\varphi}{\partial y^{m}\partial y^{j}}+\frac{\partial a_{0i}}{\partial y^{m}}\frac{\partial^{2}\varphi}{\partial y^{k}\partial y^{i}}=F(\mathbf{a}_{0}).

Next use the equation

\left(\frac{\partial^{\Theta}}{\partial y^{\Theta}}\mathcal{I}_{2}\right)_{\alpha}=0

for $|\Theta|=1$ and $\alpha=2,3$ and substitute $a_{11}$ and $\partial^{2}a_{01}$ using (13) and (14), we end up with a transport equation on $\vartheta$

(15)

\left(\frac{\partial}{\partial\tau}+A\right)\left(\begin{array}[]{c}\frac{\partial a_{02}}{\partial y^{1}}\\ \frac{\partial a_{02}}{\partial y^{2}}\\ \frac{\partial a_{02}}{\partial y^{3}}\\ \frac{\partial a_{03}}{\partial y^{1}}\\ \frac{\partial a_{03}}{\partial y^{2}}\\ \frac{\partial a_{03}}{\partial y^{3}}\end{array}\right)=F(\mathbf{a}_{0},\partial a_{01}).

Here $A$ is a $6\times 6$ matrix depending on $\lambda,\mu,\rho$ and $\varphi$ . Solving the above transport equation, one can determine $\frac{\partial^{\Theta}a_{0\alpha}}{\partial y^{\Theta}}|_{\vartheta}$ for any $|\Theta|=1$ and $\alpha=2,3$ . Using (13) and (14) again, we can determine $\partial^{2}a_{01}|_{\vartheta}$ and $a_{11}|_{\vartheta}$ .

Now we have already determined $a_{01},\partial a_{01},\partial^{2}a_{01},a_{0\alpha},\partial a_{0\alpha},a_{11}$ on $\vartheta$ , $\alpha=2,3$ . Then we can use the equations

(\mathcal{I}_{3})_{i=2,3}=0,\quad(\partial\mathcal{I}_{2})_{i=1}=0,\quad\partial^{3}\mathcal{I}_{1}=0,\quad(\partial^{2}\mathcal{I}_{2})_{i=2,3}=0

to determine $\partial^{3}a_{01},\partial^{2}a_{0\alpha},\partial a_{11},a_{1\alpha}$ for $\alpha=2,3$ .

Continuing with the process, we can have (7) satisfied, and finish our construction of Gaussian beam solutions for S-waves.

2.4. P-waves

For the construction of P-waves, we take $\varphi$ such that

\frac{\partial^{\Theta}}{\partial y^{\Theta}}(\mathcal{S}\varphi)=0\text{ on }\vartheta,

for any $|\Theta|\leq N$ , where $\mathcal{S}\varphi=(\lambda+2\mu)|\nabla\varphi|^{2}-\rho(\partial_{t}\varphi)^{2}$ . The phase function $\varphi$ can be constructed similarly as for the S-waves. Now we denote $g=g_{S}$ .

Take

\mathbf{a}_{0}=A_{P}\nabla\varphi.

Then the equation (7) is satisfied for $k=1$ since $\mathcal{I}_{1}=(\mathcal{S}\varphi)A_{P}\nabla\varphi$ .

First consider the equation $(\mathcal{I}_{2})_{i}=0$ for $i=1$ . We obtain

\begin{split}\frac{1}{\sqrt{2}}\left(\rho\partial^{2}_{t}\varphi-c_{P}^{-2}(\lambda+2\mu)\partial_{s}^{2}\varphi-c_{P}^{-2}(\lambda+2\mu)\sum_{\alpha=2}^{3}\frac{\partial^{2}\varphi}{\partial y^{\alpha}\partial y^{\alpha}}\right)A_{P}&\\ -\left(\rho\partial_{t}A_{P}+c_{P}^{-2}(\lambda+2\mu)\partial_{s}A_{P}\right)&\\ +\frac{1}{2}(\lambda+2\mu)c_{P}^{-3}\frac{\partial c_{P}}{\partial s}A_{P}-\frac{1}{2}c_{P}^{-2}\partial_{s}(\lambda+2\mu)A_{P}+\mathrm{i}\frac{1}{2}b_{1}[\rho-c_{P}^{-2}(\lambda+2\mu)]&\\ -\sqrt{2}\left(c_{P}^{-2}(\lambda+2\mu)\frac{\partial^{2}\varphi}{\partial s^{2}}+\rho\frac{\partial^{2}\varphi}{\partial s\partial t}\right)A_{P}&=0.\end{split}

We have proved in [41] that

c_{P}^{-2}(\lambda+2\mu)\frac{\partial^{2}\varphi}{\partial s^{2}}+\rho\frac{\partial^{2}\varphi}{\partial s\partial t}=\sqrt{2}\rho\frac{\partial}{\partial\tau}(\partial_{s}\varphi)=0,\quad\text{on }\vartheta,

since $\partial_{s}\varphi=\frac{1}{\sqrt{2}}$ is constant along $\vartheta$ . Using the fact that $\rho-c_{P}^{-2}(\lambda+2\mu)=0$ , the above equation can be rewritten as a transport equation

\mathcal{T}A_{P}=0

where

\mathcal{T}=2\frac{\partial}{\partial\tau}+\left[\frac{1}{\lambda+2\mu}\frac{\partial(\lambda+2\mu)}{\partial\tau}-c_{P}^{-1}\frac{\partial c_{P}}{\partial\tau}+\frac{1}{\det(Y_{P})}\frac{\partial\det(Y_{P})}{\partial\tau}\right].

Then $A_{P}|_{\vartheta}$ can be determined by solving this transport equation. In particular, we can take

A_{P}(\tau)=c\det(Y_{P}(\tau))^{-1/2}c_{P}(\tau,0)^{-1/2}\rho(\tau,0)^{-1/2},

where the constant $c$ depends on the initial value of $A_{P}$ .

Using the equation $(\mathcal{I}_{2})_{i}=0$ with $i=2,3$ , we end up with

(16)

\mathrm{i}(\rho-c_{P}^{-2}\mu)a_{1\alpha}-c_{P}^{-2}(\lambda+\mu)\frac{\partial A_{P}}{\partial y^{\alpha}}=F(A_{P}),

for $\alpha=2,3$ . Notice that $\rho-c_{P}^{-2}\mu\neq 0$ . Next use the equation

\left(\frac{\partial^{\Theta}}{\partial y^{\Theta}}\mathcal{I}_{2}\right)_{\alpha}=0

for $|\Theta|=1$ and $\alpha=1$ and substitute $a_{1\alpha}$ , $\alpha=2,3$ using (16), we end up with a transport equation

(17)

\left(\frac{\partial}{\partial\tau}+A\right)\left(\begin{array}[]{c}\frac{\partial A_{P}}{\partial y^{1}}\\ \frac{\partial A_{P}}{\partial y^{2}}\\ \frac{\partial A_{P}}{\partial y^{3}}\end{array}\right)=F(A_{P}).

Solving the above transport equation, one can determine $\frac{\partial^{\Theta}A_{P}}{\partial y^{\Theta}}|_{\vartheta}$ for any $|\Theta|=1$ . Using (16) again, we can determine $a_{1\alpha}|_{\vartheta}$ for $\alpha=2,3$ .

Similar as the construction for S-waves above, we can finish the construction of the P-wave Gaussian beam solutions such that (7) is satisfied.

3. Gaussian beams with reflections

We will discuss the reflection of Gaussian beams at the boundary in this section. The reflection condition we considered here is the traction-free boundary condition, which is natural in practice. The reflection of (real) geometric optics solutions is analyzed in [37]. We only give detailed characterization for the case where the boundary is locally flat, and non-flat case will not require much more work by flattening it [37]. Assume near a fixed point $x_{0}$ , $\partial\Omega$ is locally expressed as $x_{3}=0$ and $\Omega$ is represented as $\{x_{3}<0\}$ . We are looking for $u_{\varrho}$ as a sum of two solutions $u_{\varrho}=u_{\varrho}^{+}+u_{\varrho}^{-}$ , where $u_{\varrho}^{+}$ is the incident wave and $u_{\varrho}^{-}$ is the reflected wave generated by $u_{\varrho}^{+}$ hitting the boundary. We remark here that the reflected wave $u_{\varrho}^{-}$ can undergo further reflections, and after multiple reflections the behavior of the waves would become very complicated, but we only need to discuss single reflection.

3.1. P- incident waves

We first consider the case when the incident wave is a P-wave. The reflected wave is a combination of P- and S- waves. Assume the incident wave

u^{+}_{\varrho}=\sum_{k=0}^{N+1}\varrho^{-k}\mathbf{a}_{P,k}^{+}\chi^{+}e^{\mathrm{i}\varrho\varphi_{P}^{+}}

is a Gaussian beam traveling along the null-geodesic $\vartheta^{+}$ as constructed as in Section 2.4, and $\vartheta^{+}$ intersects with the boundary at the point $p\in(0,T)\times\partial\Omega$ . We seek for reflected waves as

(18)

u^{-}_{\varrho}=u^{-}_{P,\varrho}+u^{-}_{S,\varrho}:=\sum_{k=0}^{N+1}\varrho^{-k}\mathbf{a}_{P,k}^{-}\chi^{-}_{P}e^{\mathrm{i}\varrho\varphi_{P}^{-}}+\sum_{k=0}^{N+1}\varrho^{-k}\mathbf{a}_{S,k}^{-}\chi^{-}_{S}e^{\mathrm{i}\varrho\varphi_{S}^{-}},

where $u^{-}_{P,\varrho},u^{-}_{S,\varrho}$ are also Gaussian beam solutions representing P- and S- waves respectively. Here $\mathbf{a}^{\pm}_{\bullet,k}=(a^{\pm}_{\bullet,k1},a^{\pm}_{\bullet,k2},a^{\pm}_{\bullet,k3})$ . The reflected wave $u^{-}_{\varrho}$ is generated such that $u_{\varrho}=u^{+}_{\varrho}+u^{-}_{\varrho}$ satisfies the Neumann boundary condition

\mathcal{N}_{\lambda,\mu}u_{\varrho}:=\nu\cdot S^{L}(x,u_{\varrho})|_{\mathbb{R}\times\partial\Omega}=0,

up to order $N$ , at the point $p$ .

Assume $\nabla\varphi_{P}^{+}(p)=\xi^{+}_{P}=(\xi^{+}_{P,1},\xi^{+}_{P,2},\xi^{+}_{P,3})$ , $\nabla\varphi_{P}^{-}(p)=\xi^{-}_{P}=(\xi^{-}_{P,1},\xi^{-}_{P,2},\xi^{-}_{P,3})$ , $\nabla\varphi_{S}^{-}(p)=\xi^{-}_{S}=(\xi^{-}_{S,1},\xi^{-}_{S,2},\xi^{-}_{S,3})$ , where $|\xi^{+}_{P}|=|\xi^{-}_{P}|=\frac{1}{c_{P}(p)}$ and $|\xi^{-}_{S}|=\frac{1}{c_{S}(p)}$ , $\xi^{+}_{P,3}>0$ . By Snell’s law, we have

\xi^{+}_{P,1}=\xi^{-}_{P,1}=\xi^{-}_{S,1},\quad\xi^{+}_{P,2}=\xi^{-}_{P,2}=\xi^{-}_{S,2},\quad\xi^{-}_{P,3}=-\xi^{-}_{P,3},\quad\xi^{-}_{S,3}=-\sqrt{c_{S}^{-2}-(\xi^{+}_{P,1})^{2}-(\xi^{+}_{P,2})^{2}}.

To simplify the notation, we denote

\xi^{+}_{P}=(\xi_{1},\xi_{2},\xi_{3}),\quad\xi^{-}_{P}=(\xi_{1},\xi_{2},-\xi_{3}),\quad\xi^{-}_{S}=(\xi_{1},\xi_{2},-\xi_{S,3}).

The phase functions $\varphi^{-}_{P}$ and $\varphi^{-}_{S}$ can be constructed as in previous section such that $\varphi^{-}_{P}|_{\mathbb{R}\times\partial\Omega}=\varphi^{-}_{S}|_{\mathbb{R}\times\partial\Omega}=\varphi^{+}_{P}|_{\mathbb{R}\times\partial\Omega}$ at $p$ up to order $N$ . Then note that we have the following asymptotics

\mathcal{N}_{\lambda,\mu}u_{\varrho}=e^{\mathrm{i}\varrho\varphi}(\varrho\mathcal{N}_{0}+\mathcal{N}_{1}+\varrho^{-1}\mathcal{N}_{2}+\cdots),

where $\varphi=\varphi^{+}_{P}|_{\mathbb{R}\times\partial\Omega}$ . We need

(19)

\partial^{\Theta}\mathcal{N}_{k}=0\text{ at }p,\quad\text{for }k\geq 0,\,|\Theta|\leq N.

By calculation, we obtain

\begin{split}\mathcal{N}_{0}=&\left(\begin{array}[]{ccc}\mu\partial_{3}\varphi_{P}^{+}&0&\mu\partial_{1}\varphi_{P}^{+}\\ 0&\mu\partial_{3}\varphi_{P}^{+}&\mu\partial_{2}\varphi_{P}^{+}\\ \lambda\partial_{1}\varphi_{P}^{+}&\lambda\partial_{2}\varphi_{P}^{+}&(\lambda+2\mu)\partial_{3}\varphi_{P}^{+}\end{array}\right)\left(\begin{array}[]{c}a_{P,01}^{+}\\ a_{P,02}^{+}\\ a_{P,03}^{+}\end{array}\right)\\ &+\left(\begin{array}[]{ccc}\mu\partial_{3}\varphi_{P}^{-}&0&\mu\partial_{1}\varphi_{P}^{-}\\ 0&\mu\partial_{3}\varphi_{P}^{-}&\mu\partial_{2}\varphi_{P}^{-}\\ \lambda\partial_{1}\varphi_{P}^{-}&\lambda\partial_{2}\varphi_{P}^{-}&(\lambda+2\mu)\partial_{3}\varphi_{P}^{-}\end{array}\right)\left(\begin{array}[]{c}a_{P,01}^{-}\\ a_{P,02}^{-}\\ a_{P,03}^{-}\end{array}\right)\\ &+\left(\begin{array}[]{ccc}\mu\partial_{3}\varphi_{S}^{-}&0&\mu\partial_{1}\varphi_{S}^{-}\\ 0&\mu\partial_{3}\varphi_{S}^{-}&\mu\partial_{2}\varphi_{S}^{-}\\ \lambda\partial_{1}\varphi_{S}^{-}&\lambda\partial_{2}\varphi_{S}^{-}&(\lambda+2\mu)\partial_{3}\varphi_{S}^{-}\end{array}\right)\left(\begin{array}[]{c}a_{S,01}^{-}\\ a_{S,02}^{-}\\ a_{S,03}^{-}\end{array}\right).\end{split}

Using the fact

\mathbf{a}_{P,0}^{-}=A_{P}^{-}\nabla\varphi_{P}^{-},\quad\mathbf{a}_{P,0}^{+}=A_{P}^{+}\nabla\varphi_{P}^{+},

we have

(20)

\begin{split}\mathcal{N}_{0}=&\left(\begin{array}[]{ccc}\mu\partial_{3}\varphi_{P}^{-}&0&\mu\partial_{1}\varphi_{P}^{-}\\ 0&\mu\partial_{3}\varphi_{P}^{-}&\mu\partial_{2}\varphi_{P}^{-}\\ \lambda\partial_{1}\varphi_{P}^{-}&\lambda\partial_{2}\varphi_{P}^{-}&(\lambda+2\mu)\partial_{3}\varphi_{P}^{-}\end{array}\right)\left(\begin{array}[]{c}\partial_{1}\varphi_{P}^{-}\\ \partial_{2}\varphi_{P}^{-}\\ \partial_{3}\varphi_{P}^{-}\end{array}\right)A_{P}^{-}\\ &+\left(\begin{array}[]{ccc}\mu\partial_{3}\varphi_{S}^{-}&0&\mu\partial_{1}\varphi_{S}^{-}\\ 0&\mu\partial_{3}\varphi_{S}^{-}&\mu\partial_{2}\varphi_{S}^{-}\\ \lambda\partial_{1}\varphi_{S}^{-}&\lambda\partial_{2}\varphi_{S}^{-}&(\lambda+2\mu)\partial_{3}\varphi_{S}^{-}\end{array}\right)\left(\begin{array}[]{c}a_{S,01}^{-}\\ a_{S,02}^{-}\\ a_{S,03}^{-}\end{array}\right)\\ &+\left(\begin{array}[]{ccc}\mu\partial_{3}\varphi_{P}^{+}&0&\mu\partial_{1}\varphi_{P}^{+}\\ 0&\mu\partial_{3}\varphi_{P}^{+}&\mu\partial_{2}\varphi_{P}^{+}\\ \lambda\partial_{1}\varphi_{P}^{+}&\lambda\partial_{2}\varphi_{P}^{+}&(\lambda+2\mu)\partial_{3}\varphi_{P}^{+}\end{array}\right)\left(\begin{array}[]{c}\partial_{1}\varphi_{P}^{+}\\ \partial_{2}\varphi_{P}^{+}\\ \partial_{3}\varphi_{P}^{+}\end{array}\right)A_{P}^{+}.\end{split}

Considering also the identity

(21)

\langle\mathbf{a}_{S,0}^{-},\nabla\varphi_{S}^{-}\rangle=0,

we have the following equations at point $p$ :

\begin{split}\left(\begin{array}[]{ccc}-\mu\xi_{3}&0&\mu\xi_{1}\\ 0&-\mu\xi_{3}&\mu\xi_{2}\\ \lambda\xi_{1}&\lambda\xi_{2}&-(\lambda+2\mu)\xi_{3}\end{array}\right)\left(\begin{array}[]{c}\xi_{1}\\ \xi_{2}\\ -\xi_{3}\end{array}\right)A_{P}^{-}(p)\\ +\left(\begin{array}[]{ccc}-\mu\xi_{S,3}&0&\mu\xi_{1}\\ 0&-\mu\xi_{S,3}&\mu\xi_{2}\\ \lambda\xi_{1}&\lambda\xi_{2}&-(\lambda+2\mu)\xi_{S,3}\end{array}\right)\left(\begin{array}[]{c}a_{S,01}^{-}(p)\\ a_{S,02}^{-}(p)\\ a_{S,03}^{-}(p)\end{array}\right)\\ =-\left(\begin{array}[]{ccc}\mu\xi_{3}&0&\mu\xi_{1}\\ 0&\mu\xi_{3}&\mu\xi_{2}\\ \lambda\xi_{1}&\lambda\xi_{2}&(\lambda+2\mu)\xi_{3}\end{array}\right)\left(\begin{array}[]{c}\xi_{1}\\ \xi_{2}\\ \xi_{3}\end{array}\right)A_{P}^{+}(p).\end{split}

and

(\xi_{1},\xi_{2},-\xi_{S,3})\left(\begin{array}[]{c}a_{S,01}^{-}(p)\\ a_{S,02}^{-}(p)\\ a_{S,03}^{-}(p)\end{array}\right)=0.

The equations above can be formulated into the following linear system for $(A_{P}^{-},a_{S,01}^{-},a_{S,02}^{-},a_{S,03}^{-})$ at $p$ :

(22)

\begin{split}M_{P}(\xi)\left(\begin{array}[]{c}A_{P}^{-}(p)\\ a_{S,01}^{-}(p)\\ a_{S,02}^{-}(p)\\ a_{S,03}^{-}(p)\end{array}\right):=&\left(\begin{array}[]{cccc}-2\mu\xi_{1}\xi_{3}&-\mu\xi_{S,3}&0&\mu\xi_{1}\\ -2\mu\xi_{2}\xi_{3}&0&-\mu\xi_{S,3}&\mu\xi_{2}\\ \rho-2\mu(\xi_{1}^{2}+\xi_{2}^{2})&\lambda\xi_{1}&\lambda\xi_{2}&-(\lambda+2\mu)\xi_{S,3}\\ 0&\xi_{1}&\xi_{2}&-\xi_{S,3}\end{array}\right)\left(\begin{array}[]{c}A_{P}^{-}(p)\\ a_{S,01}^{-}(p)\\ a_{S,02}^{-}(p)\\ a_{S,03}^{-}(p)\end{array}\right)\\ =&-\left(\begin{array}[]{c}-2\mu\xi_{1}\xi_{3}A_{P}^{+}(p)\\ -2\mu\xi_{2}\xi_{3}A_{P}^{+}(p)\\ \rho-2\mu(\xi_{1}^{2}+\xi_{2}^{2})A_{P}^{+}(p)\\ 0\end{array}\right).\end{split}

Lemma 1.

The matrix $M_{P}(\xi)$ is invertible.

Proof.

To see that the matrix $M_{P}(\xi)$ is invertible, without loss of generality, we assume $\xi_{2}=0$ . Consider the homogeneous equation

(23)

M(\xi)\left(\begin{array}[]{c}A_{P}\\ a_{S,1}\\ a_{S,2}\\ a_{S,3}\end{array}\right)=0.

By the last equation, we can assume that

\left(\begin{array}[]{c}a_{S,1}\\ a_{S,2}\\ a_{S,3}\end{array}\right)=A_{SH}\mathbf{e}_{SH}+A_{SV}\mathbf{e}_{SV}

where

\mathbf{e}_{SH}=\left(\begin{array}[]{c}0\\ -c_{S}^{-1}\\ 0\end{array}\right)

\mathbf{e}_{SV}=\left(\begin{array}[]{c}\xi_{S,3}\\ 0\\ \xi_{1}\end{array}\right).

The homogeneous equation simplifies into

M^{\prime}_{P}(\xi)(A_{P},A_{SV},A_{SH})^{T}=0,

where

M_{P}^{\prime}(\xi)=\left(\begin{array}[]{ccc}-2\mu\xi_{1}\xi_{3}&-\rho+2\mu\xi_{1}^{2}&0\\ 0&0&-\mu c_{S}^{-1}\xi_{S,3}\\ \rho-2\mu\xi_{1}^{2}&-2\mu\xi_{1}\xi_{S,3}&0\end{array}\right).

One immediately obtain that $A_{SH}=0$ . Then we only need to consider the matrix

\left(\begin{array}[]{cc}-2\mu\xi_{1}\xi_{3}&\mu(\xi_{1}^{2}-\xi_{S,3}^{2})\\ -\mu(\xi_{1}^{2}-\xi_{S,3}^{2})&-2\mu\xi_{1}\xi_{S,3}\end{array}\right),

whose determinant is

\mu^{2}(4\xi_{1}^{2}\xi_{3}\xi_{S,3}+\xi_{1}^{4}-2\xi_{1}^{2}\xi_{S,3}^{2}+\xi_{S,3}^{4})>0.

This is because that $\xi_{1}^{2}\xi_{3}\xi_{S,3}<\xi_{1}^{2}\xi_{S,3}^{2}$ , which in turn yields

0<4\xi_{1}^{2}\xi_{3}\xi_{S,3}-2\xi_{1}^{2}\xi_{S,3}^{2}<2\xi_{1}^{2}\xi_{S,3}^{2}\leq\xi_{1}^{4}+\xi_{S,3}^{4}.

This shows that the equation (23) has has only zero solutions, and therefore $M(\xi)$ is invertible. ∎

By the above lemma, we can solve the linear system (22) to determine $(A_{P}^{-},a_{S,01}^{-},a_{S,02}^{-},a_{S,03}^{-})$ at the point $p$ .

Next, we determine the tangential derivatives of $(A_{P}^{-},a_{S,01}^{-},a_{S,02}^{-},a_{S,03}^{-})$ at $p$ . Use the fact that $\mathcal{N}_{0}=0$ needs to be satisfied up to order $N$ at $p$ , that is $\partial_{t,x}^{\alpha}\mathcal{N}_{1}(p)=0$ , for $\alpha=(\alpha_{0},\alpha_{1},\alpha_{2},0)$ , $|\alpha|\leq N$ . Let us first consider, for example, $(\partial_{1}A_{P}^{-},\partial_{1}a_{S,01}^{-},\partial_{1}a_{S,02}^{-},\partial_{1}a_{S,03}^{-})$ . Taking the first order derivative in $x_{1}$ of (20) and (21), we obtain a linear system

M(\xi)\left(\begin{array}[]{c}\partial_{1}A_{P}^{-}(p)\\ \partial_{1}a_{S,01}^{-}(p)\\ \partial_{1}a_{S,02}^{-}(p)\\ \partial_{1}a_{S,03}^{-}(p)\end{array}\right)=F(A_{P}^{-}(p),\mathbf{a}_{S,0}^{-}(p)).

Here and below we suppress the dependence of $F$ on $\mathbf{a}^{+}_{P,k}$ and $\varphi^{+}_{P},\varphi^{-}_{P},\varphi^{-}_{S}$ (as well as their derivates). Because the matrix $M(\xi)$ is invertible, we can determine $(\partial_{1}A_{P}^{-},\partial_{1}a_{S,01}^{-},\partial_{1}a_{S,02}^{-},\partial_{1}a_{S,03}^{-})$ at the point $p$ . Continue with this process, one can determine $(\partial^{\alpha}A_{P}^{-},\partial^{\alpha}a_{S,01}^{-},\partial^{\alpha}a_{S,02}^{-},\partial^{\alpha}a_{S,03}^{-})(p)$ for any $\alpha=(\alpha_{0},\alpha_{1},\alpha_{2},0)$ , $|\alpha|\leq N$ .

Next, consider the determination of $\mathbf{a}_{P,1}^{-}$ and $\mathbf{a}_{S,1}^{-}$ at the point $p$ . Notice that

(24)

\begin{split}\mathcal{N}_{1}=&\left(\begin{array}[]{ccc}\mu\partial_{3}\varphi_{P}^{-}&0&\mu\partial_{1}\varphi_{P}^{-}\\ 0&\mu\partial_{3}\varphi_{P}^{-}&\mu\partial_{2}\varphi_{P}^{-}\\ \lambda\partial_{1}\varphi_{P}^{-}&\lambda\partial_{2}\varphi_{P}^{-}&(\lambda+2\mu)\partial_{3}\varphi_{P}^{-}\end{array}\right)\left(\begin{array}[]{c}a_{P,11}^{-}\\ a_{P,12}^{-}\\ a_{P,13}^{-}\end{array}\right)\\ &+\left(\begin{array}[]{ccc}\mu\partial_{3}\varphi_{S}^{-}&0&\mu\partial_{1}\varphi_{S}^{-}\\ 0&\mu\partial_{3}\varphi_{S}^{-}&\mu\partial_{2}\varphi_{S}^{-}\\ \lambda\partial_{1}\varphi_{S}^{-}&\lambda\partial_{2}\varphi_{S}^{-}&(\lambda+2\mu)\partial_{3}\varphi_{S}^{-}\end{array}\right)\left(\begin{array}[]{c}a_{S,11}^{-}\\ a_{S,12}^{-}\\ a_{S,13}^{-}\end{array}\right)\\ &+F(\mathbf{a}^{-}_{P,0},\mathbf{a}^{-}_{S,0},\partial\mathbf{a}^{-}_{P,0},\partial\mathbf{a}^{-}_{S,0})\end{split}

Note that the full first-order derivatives of $\mathbf{a}^{-}_{P,0},\mathbf{a}^{-}_{S,0}$ at point $p$ , $\partial\mathbf{a}^{-}_{P,0}(p),\partial\mathbf{a}^{-}_{S,0}(p)$ (not only the derivatives in $(t,x_{1},x_{2})$ ), are determined since the transport equations for $A_{P}^{-}$ and $\mathbf{a}^{-}_{S,0}$ are satisfied. The equation $\mathcal{N}_{1}(p)=0$ can be rewritten as

(25)

\begin{split}\left(\begin{array}[]{ccc}-\mu\xi_{3}&0&\mu\xi_{1}\\ 0&-\mu\xi_{3}&\mu\xi_{2}\\ \lambda\xi_{1}&\lambda\xi_{2}&-(\lambda+2\mu)\xi_{3}\end{array}\right)\left(\begin{array}[]{c}a_{P,11}^{-}(p)\\ a_{P,12}^{-}(p)\\ a_{P,13}^{-}(p)\end{array}\right)\\ +\left(\begin{array}[]{ccc}-\mu\xi_{S,3}&0&\mu\xi_{1}\\ 0&-\mu\xi_{S,3}&\mu\xi_{2}\\ \lambda\xi_{1}&\lambda\xi_{2}&-(\lambda+2\mu)\xi_{S,3}\end{array}\right)\left(\begin{array}[]{c}a_{S,11}^{-}(p)\\ a_{S,12}^{-}(p)\\ a_{S,13}^{-}(p)\end{array}\right)=F,\end{split}

where $F$ have already been computed in the previous steps. Now, we have $3$ equations for the $6$ unknowns. The extra $3$ equations come from equations (13) and (16), which can be rewritten as

(26)

(\xi_{1},\xi_{2},-\xi_{S,3})\cdot(a_{S,11}^{-}(p),a_{S,12}^{-}(p),a_{S,13}^{-}(p))=F_{1}(\mathbf{a}^{-}_{P,0,}(p),\partial\mathbf{a}^{-}_{P,0,}(p),\mathbf{a}^{-}_{S,0}(p),\partial\mathbf{a}^{-}_{S,0,}(p)),

and

(27)

\begin{split}(\xi_{2},-\xi_{1},0)\cdot(a_{P,11}^{-}(p),a_{P,12}^{-}(p),a_{P,13}^{-})(p)=F_{2}(\mathbf{a}^{-}_{P,0}(p),\partial\mathbf{a}^{-}_{P,0}(p),\mathbf{a}^{-}_{S,0}(p),\partial\mathbf{a}^{-}_{S,0}(p)),\\ (\xi_{1}\xi_{3},\xi_{2}\xi_{3},\xi_{1}^{2}+\xi_{2}^{2})\cdot(a_{P,11}^{-}(p),a_{P,12}^{-}(p),a_{P,13}^{-}(p))=F_{3}(\mathbf{a}^{-}_{P,0}(p),\partial\mathbf{a}^{-}_{P,0}(p),\mathbf{a}^{-}_{S,0}(p),\partial\mathbf{a}^{-}_{S,0}(p)).\end{split}

Using (27), we can reduce the linear system (25) to

\begin{split}\left(\begin{array}[]{ccc}-\mu\xi_{3}&0&\mu\xi_{1}\\ 0&-\mu\xi_{3}&\mu\xi_{2}\\ \lambda\xi_{1}&\lambda\xi_{2}&-(\lambda+2\mu)\xi_{3}\end{array}\right)\left(\begin{array}[]{c}\xi_{1}\\ \xi_{2}\\ -\xi_{3}\end{array}\right)A_{P,1}^{-}(p)\\ +\left(\begin{array}[]{ccc}-\mu\xi_{S,3}&0&\mu\xi_{1}\\ 0&-\mu\xi_{S,3}&\mu\xi_{2}\\ \lambda\xi_{1}&\lambda\xi_{2}&-(\lambda+2\mu)\xi_{S,3}\end{array}\right)\left(\begin{array}[]{c}a_{S,11}^{-}(p)\\ a_{S,12}^{-}(p)\\ a_{S,13}^{-}(p)\end{array}\right)=F,\end{split}

where $A_{P,1}^{-}(p)=(a_{P,11}^{-}(p),a_{P,12}^{-}(p),a_{P,13}^{-}(p))\cdot(\xi_{1},\xi_{2},-\xi_{3})$ . Together with (26), we now have a system

M_{P}(\xi)(A_{P,1}^{-}(p),a_{S,11}^{-}(p),a_{S,12}^{-}(p),a_{S,13}^{-}(p))^{T}=F.

By Lemma 1, the matrix $M_{P}(\xi)$ is invertible. Then, $a_{P,11}^{-},a_{P,12}^{-},a_{P,13}^{-},a_{S,11}^{-},a_{S,12}^{-},a_{S,13}^{-}$ at point $p$ can be determined.

To determine the first order derivative (for example in $x_{1}$ ) of $\mathbf{a}_{P,1}^{-}$ and $\mathbf{a}_{S,1}^{-}$ at $p$ , we take the $\partial_{x_{1}}$ derivative of (25), and also use the equations

\begin{split}(\xi_{1},\xi_{2},-\xi_{S,3})\cdot(\partial_{x_{1}}a_{S,11}^{-}(p),\partial_{x_{1}}a_{S,12}^{-}(p),\partial_{x_{1}}a_{S,13}^{-}(p))=F_{1},\\ (\xi_{2},-\xi_{1},0)\cdot(\partial_{x_{1}}a_{P,11}^{-}(p),\partial_{x_{1}}a_{P,12}^{-}(p),\partial_{x_{1}}a_{P,13}^{-}(p))=F_{2},\\ (\xi_{1}\xi_{3},\xi_{2}\xi_{3},\xi_{1}^{2}+\xi_{2}^{2})\cdot(\partial_{x_{1}}a_{P,11}^{-}(p),\partial_{x_{1}}a_{P,12}^{-}(p),\partial_{x_{1}}a_{P,13}^{-}(p))=F_{3},\end{split}

which come from $\partial_{x_{1}}(\mathcal{I}_{2})_{i}(p)=0$ (more precisely, $i=2,3$ for P-wave, and $i=1$ for S-wave). Consequently, we end up with a linear system as

\left(\begin{array}[]{cccccc}-\mu\xi_{3}&0&\mu\xi_{1}&-\mu\xi_{S,3}&0&\mu\xi_{1}\\ 0&-\mu\xi_{3}&\mu\xi_{2}&0&-\mu\xi_{S,3}&\mu\xi_{2}\\ \lambda\xi_{1}&\lambda\xi_{2}&-(\lambda+2\mu)\xi_{3}&\lambda\xi_{1}&\lambda\xi_{2}&-(\lambda+2\mu)\xi_{S,3}\\ 0&0&0&\xi_{1}&\xi_{2}&-\xi_{S,3}\\ \xi_{2}&-\xi_{1}&0&0&0&0\\ \xi_{1}\xi_{3}&\xi_{2}\xi_{3}&\xi_{1}^{2}+\xi_{2}^{2}&0&0&0\end{array}\right)\left(\begin{array}[]{c}\partial_{1}a_{P,11}^{-}(p)\\ \partial_{1}a_{P,12}^{-}(p)\\ \partial_{1}a_{P,13}^{-}(p)\\ \partial_{1}a_{S,11}^{-}(p)\\ \partial_{1}a_{S,12}^{-}(p)\\ \partial_{1}a_{S,13}^{-}(p)\end{array}\right)=F.

By similar consideration as above, this linear system is solvable, and thus we can determine $\partial_{1}a_{P,11}^{-},\partial_{1}a_{P,12}^{-},\partial_{1}a_{P,13}^{-},\partial_{1}a_{S,11}^{-},\partial_{1}a_{S,12}^{-},\partial_{1}a_{S,13}^{-}$ at $p$ .

Continuing with this process, we can determine $\partial^{\alpha}\mathbf{a}^{-}_{P,k}(p)$ and $\partial^{\alpha}\mathbf{a}^{-}_{S,k}(p)$ for any $k\leq N$ and $|\alpha|\leq N$ . Now we finish the construction of the reflected waves.

3.2. S- incident waves

In this section, we assume that the incident wave is S-wave, i.e.,

u^{+}_{\varrho}=e^{\mathrm{i}\varrho\varphi_{S}^{+}}\chi\left(\sum_{k=0}^{N+1}\varrho^{-k}\mathbf{a}_{S,k}^{+}\right).

We construct reflected waves given by (18). Denote

\xi^{+}_{S}:=\xi=(\xi_{1},\xi_{2},\xi_{3}),\quad\xi^{-}_{S}=(\xi_{1},\xi_{2},-\xi_{3}),\quad\xi^{-}_{P}=(\xi_{1},\xi_{2},-\xi_{P,3}).

where

(28)

\xi_{P,3}=\sqrt{c_{P}^{-2}-\xi_{1}^{2}-\xi_{2}^{2}}.

Then the phase functions satisfy

\nabla\varphi_{S}^{+}(p)=\xi^{+}_{S},\quad\nabla\varphi_{P}^{-}(p)=\xi^{-}_{P},\quad\nabla\varphi_{S}^{-}(p)=\xi^{-}_{S}.

First we consider the important case, for which $\xi_{1}^{2}+\xi_{2}^{2}<c_{P}^{-2}$ , so that $\xi_{P,3}$ is real and there is no evanescent wave. The reflected P- and S- waves are both progressing waves. The equation $\mathcal{N}_{0}(p)=0$ would give a system

M_{S}(\xi)\left(\begin{array}[]{c}A_{P}^{-}(p)\\ a_{S,01}^{-}(p)\\ a_{S,02}^{-}(p)\\ a_{S,03}^{-}(p)\end{array}\right)=F,

where

M_{S}(\xi)=\left(\begin{array}[]{cccc}-2\mu\xi_{1}\xi_{P,3}&-\mu\xi_{3}&0&\mu\xi_{1}\\ -2\mu\xi_{2}\xi_{P,3}&0&-\mu\xi_{3}&\mu\xi_{2}\\ \rho-2\mu(\xi_{1}^{2}+\xi_{2}^{2})&\lambda\xi_{1}&\lambda\xi_{2}&-(\lambda+2\mu)\xi_{3}\\ 0&\xi_{1}&\xi_{2}&-\xi_{3}\end{array}\right).

The matrix $M_{S}(\xi)=M_{P}((\xi_{1},\xi_{2},\xi_{P,3}))$ is invertible. Therefore, similar as in previous section, we can determine $\partial^{\alpha}\mathbf{a}^{-}_{P,k}(p)$ and $\partial^{\alpha}\mathbf{a}^{-}_{S,k}(p)$ for any $k\leq N$ and $|\alpha|\leq N$ .

Evanescent waves. Now we consider the case $\xi_{1}^{2}+\xi_{2}^{2}>c_{P}^{-2}.$ In this case, $\xi_{P,3}$ is taken by

\xi_{P,3}=\mathrm{i}\sqrt{\xi_{1}^{2}+\xi_{2}^{2}-c_{P}^{-2}}.

To be consistent with the formula (28), we can choose the square root function $\sqrt{\,\cdot\,}$ such that $\Im\sqrt{z}>0$ for $z\in\mathbb{C}\setminus[0,+\infty)$ . Then $\Re(-\mathrm{i}\xi_{P,3})>0$ , and we can then construct $\varphi_{P}^{-}$ in a neighborhood of $p$ such that $\partial_{x_{3}}(\varphi_{P}^{-})=-\xi_{P,3}$ . Then $\Re(\mathrm{i}\varphi_{P}^{-})<0$ for $x_{3}<0$ and consequently

|e^{\mathrm{i}\varrho\varphi_{P}^{-}}|=\mathcal{O}(|x_{3}|^{\infty})

near $p$ in $\mathbb{R}\times\Omega$ . Then we can construct $u_{\varrho}^{-}$ of the form (18) in a neighborhood of $p$ such that

\rho\partial_{t}^{2}u_{P,\varrho}^{-}-\nabla\cdot S^{L}(u_{P,\varrho}^{-})=0,\quad\text{up to order }N\text{ at point }p,

and now $\chi^{-}_{P}$ is compactly supported in a neighborhood of $p$ . We also refer to [37] for a discussion on evanescent waves.

3.3. Full Gaussian beam asymptotic solutions

Now it is clear to see how to construct Gaussian beam solutions to the linear elastic wave equation incorporating all reflections at the boundary $(0,T)\times\partial\Omega$ .

Assume $\vartheta:(t^{-},t^{+})\rightarrow M$ be a unit speed null-geodesic in $(M,-\mathrm{d}t^{2}+c_{\bullet}^{-2}\mathrm{d}s^{2})$ , with endpoints $\vartheta(t^{-}),\vartheta(t^{+})\in\partial M$ . Let $R_{0}\subset\partial M$ be a neighborhood of $\vartheta(t^{-})$ . Assume that $t^{-}\in(0,T)$ . If $\vartheta$ is a forward null-geodesic then $t^{+}>t^{-}$ , if $\vartheta$ is a backward null-geodesic then $t^{+}<t^{-}$ .

Fix $k$ and $K$ , by above discussions, we can take $N$ large enough and construct asymptotic solutions $u_{\varrho}$ such that

\|\mathcal{L}_{\lambda,\mu,\rho}u_{\varrho}\|_{H^{k}((0,T)\times\Omega)}=\mathcal{O}(\varrho^{-K}),

the boundary values of $u_{\varrho}$ satisfy

\|\mathcal{N}_{\lambda,\mu}u_{\varrho}\|_{H^{k}((0,T)\times\partial\Omega\setminus R_{0})}=\mathcal{O}(\varrho^{-K}).

Actually, $u_{\varrho}=u_{\varrho}^{\mathrm{inc}}+u_{\varrho}^{\mathrm{ref}}$ , where the incident wave $u_{\varrho}^{\mathrm{inc}}=u_{\varrho}^{+}$ is compactly supported in a neighborhood of $\vartheta$ . We remark here that $u_{\varrho}^{\mathrm{inc}}$ is a Gaussian beam starting from $\vartheta(t_{-})$ , and $u_{\rho}^{\mathrm{ref}}$ is generated by the reflection of the incident wave $u_{\varrho}^{\mathrm{inc}}$ at $\vartheta(t_{+})$ and subsequent reflections.

4. Proof of the main theorem

4.1. Second order linearization of the displacement-to-traction map

We first summarize the second order linearization of the displacement-to-traction map carried out in [41]. We also refer to [21, 24, 17] for the use of higher order linearization in the study of inverse problems for nonlinear equations.

Take $\epsilon_{1},\epsilon_{2}\in\mathbb{R}$ small enough, and let $u_{\epsilon}$ be the solution to the initial boundary value problem (2) with boundary value $f=\epsilon_{1}f^{(1)}+\epsilon_{2}f^{(2)}$ , then $u_{\epsilon}$ has the asymptotic expansion

u_{\epsilon}=\epsilon_{1}u^{(1)}+\epsilon_{2}u^{(2)}+\frac{1}{2}\epsilon_{1}^{2}u^{(11)}+\frac{1}{2}\epsilon_{2}^{2}u^{(22)}+\epsilon_{1}\epsilon_{2}u^{(12)}+\text{higher order terms in }\epsilon_{1},\epsilon_{2}.

Here $u^{(1)},u^{(2)}$ are solutions to the linearized equation

(29)

\begin{split}&\rho\frac{\partial^{2}u^{(j)}}{\partial t^{2}}-\nabla\cdot S^{L}(x,u^{(j)})=0,\quad(t,x)\in(0,T)\times\Omega,\\ &u^{(j)}(t,x)=f^{(j)}(t,x),\quad(t,x)\in(0,T)\times\partial\Omega,\\ &u^{(j)}(0,x)=\frac{\partial}{\partial t}u^{(j)}(0,x)=0,\quad x\in\Omega,\end{split}

and $u^{(12)}$ is the solution to the equation

(30)

\begin{split}&\rho\frac{\partial^{2}u^{(12)}}{\partial t^{2}}-\nabla\cdot S^{L}(x,u^{(12)})=\nabla\cdot G(u^{(1)},u^{(2)}),\quad(t,x)\in(0,T)\times\Omega,\\ &u^{(12)}(t,x)=0,\quad(t,x)\in(0,T)\times\partial\Omega,\\ &u^{(12)}(0,x)=\frac{\partial}{\partial t}u^{(12)}(0,x)=0,\quad x\in\Omega,\end{split}

where $G_{ij}$ comes from the second order term of $u$ in $S(x,u)$ ,

(31)

\begin{split}G_{ij}(u^{(1)},u^{(2)})=&(\lambda+\mathscr{B})\frac{\partial u^{(1)}_{m}}{\partial x_{n}}\frac{\partial u^{(2)}_{m}}{\partial x_{n}}\delta_{ij}+2\mathscr{C}\frac{\partial u^{(1)}_{m}}{\partial x_{m}}\frac{\partial u^{(2)}_{n}}{\partial x_{n}}\delta_{ij}+\mathscr{B}\frac{\partial u_{m}^{(1)}}{\partial x_{n}}\frac{\partial u_{n}^{(2)}}{\partial x_{m}}\delta_{ij}\\ &+\mathscr{B}\left(\frac{\partial u^{(1)}_{m}}{\partial x_{m}}\frac{\partial u^{(2)}_{j}}{\partial x_{i}}+\frac{\partial u^{(2)}_{m}}{\partial x_{m}}\frac{\partial u^{(1)}_{j}}{\partial x_{i}}\right)+\frac{\mathscr{A}}{4}\left(\frac{\partial u^{(1)}_{j}}{\partial x_{m}}\frac{\partial u^{(2)}_{m}}{\partial x_{i}}+\frac{\partial u^{(2)}_{j}}{\partial x_{m}}\frac{\partial u^{(1)}_{m}}{\partial x_{i}}\right)\\ &+(\lambda+\mathscr{B})\left(\frac{\partial u^{(1)}_{m}}{\partial x_{m}}\frac{\partial u^{(2)}_{i}}{\partial x_{j}}+\frac{\partial u^{(2)}_{m}}{\partial x_{m}}\frac{\partial u^{(1)}_{i}}{\partial x_{j}}\right)\\ &+\left(\mu+\frac{\mathscr{A}}{4}\right)\Bigg{(}\frac{\partial u^{(1)}_{m}}{\partial x_{i}}\frac{\partial u^{(2)}_{m}}{\partial x_{j}}+\frac{\partial u^{(2)}_{m}}{\partial x_{i}}\frac{\partial u^{(1)}_{m}}{\partial x_{j}}+\frac{\partial u^{(1)}_{i}}{\partial x_{m}}\frac{\partial u^{(2)}_{j}}{\partial x_{m}}+\frac{\partial u^{(2)}_{i}}{\partial x_{m}}\frac{\partial u^{(1)}_{j}}{\partial x_{m}}\\ &\quad\quad+\frac{\partial u^{(1)}_{i}}{\partial x_{m}}\frac{\partial u^{(2)}_{m}}{\partial x_{j}}+\frac{\partial u^{(2)}_{i}}{\partial x_{m}}\frac{\partial u^{(1)}_{m}}{\partial x_{j}}\Bigg{)}.\end{split}

We define the linear map

\begin{split}\Lambda^{{}^{\prime}}_{D}:f^{(j)}&\mapsto\nu\cdot S^{L}(x,u^{(j)})|_{(0,T)\times\partial\Omega},\\ \end{split}

and the bilinear map

\begin{split}\Lambda^{{}^{\prime\prime}}_{D}:(f^{(1)},f^{(2)})&\mapsto\left(\nu\cdot S^{L}(u^{(12)})+\nu\cdot G(u^{(1)},u^{(2)})\right)\Big{|}_{(0,T)\times\partial\Omega}.\\ \end{split}

Here $f^{(1)},f^{(2)}$ vanish near $\{t=0\}$ . We remark here that formally

\begin{split}\Lambda^{{}^{\prime}}_{D}(f^{(j)})&=\frac{\partial}{\partial\epsilon_{j}}\Lambda(\epsilon_{j}f^{(j)})|_{\epsilon_{j}=0},\\ \Lambda^{{}^{\prime\prime}}_{D}(f^{(1)},f^{(2)})&=\frac{\partial^{2}}{\partial\epsilon_{1}\partial\epsilon_{2}}\Lambda(\epsilon_{1}f^{(1)}+\epsilon_{2}f^{(2)})|_{\epsilon_{1}=\epsilon_{2}=0}.\end{split}

Therefore, we can recover $\Lambda^{\prime}_{D},\Lambda^{\prime\prime}_{D}$ from $\Lambda$ .

Assume that $u^{(0)}$ solves the initial boundary value problem for the backward elastic wave equation

(32)

\begin{split}&\rho\frac{\partial^{2}}{\partial t^{2}}u^{(0)}-\nabla\cdot S^{L}(x,u^{(0)})=0,\quad(t,x)\in(0,T)\times\Omega,\\ &u(t,x)=f^{(0)}(t,x),\quad(t,x)\in(0,T)\times\partial\Omega,\\ &u^{(0)}(T,x)=\frac{\partial}{\partial t}u^{(0)}(T,x)=0,\quad x\in\Omega.\end{split}

Using integration by parts, we have (cf. [41])

(33)

\int_{0}^{T}\int_{\partial\Omega}\Lambda^{{}^{\prime\prime}}_{D}(f^{(1)},f^{(2)})\cdot f^{(0)}\mathrm{d}S\mathrm{d}t\\ =\int_{0}^{T}\int_{\Omega}\mathcal{G}(\nabla u^{(1)},\nabla u^{(2)},\nabla u^{(0)})\mathrm{d}x\mathrm{d}t,

where

(34)

\begin{split}&\mathcal{G}(\nabla u^{(1)},\nabla u^{(2)},\nabla u^{(0)})\\ =&(\lambda+\mathscr{B})(\nabla u^{(1)}:\nabla u^{(2)})(\nabla\cdot u^{(0)})+2\mathscr{C}(\nabla\cdot u^{(1)})(\nabla\cdot u^{(2)})(\nabla\cdot u^{(0)})\\ &+\mathscr{B}(\nabla u^{(1)}:\nabla^{T}u^{(2)})(\nabla\cdot u^{(0)})\\ &+\mathscr{B}\left((\nabla\cdot u^{(1)})(\nabla u^{(2)}:\nabla^{T}u^{(0)})+(\nabla\cdot u^{(2)})(\nabla u^{(1)}:\nabla^{T}u^{(0)})\right)\\ &+\frac{\mathscr{A}}{4}\left(\frac{\partial u^{(1)}_{j}}{\partial x_{m}}\frac{\partial u^{(2)}_{m}}{\partial x_{i}}+\frac{\partial u^{(2)}_{j}}{\partial x_{m}}\frac{\partial u^{(1)}_{m}}{\partial x_{i}}\right)\frac{\partial u^{(0)}_{i}}{\partial x_{j}}\\ &+(\lambda+\mathscr{B})\left((\nabla\cdot u^{(1)})(\nabla u^{(2)}:\nabla u^{(0)})+(\nabla\cdot u^{(2)})(\nabla u^{(1)}:\nabla u^{(0)})\right)\\ &+\left(\mu+\frac{\mathscr{A}}{4}\right)\Bigg{(}\frac{\partial u^{(1)}_{m}}{\partial x_{i}}\frac{\partial u^{(2)}_{m}}{\partial x_{j}}+\frac{\partial u^{(2)}_{m}}{\partial x_{i}}\frac{\partial u^{(1)}_{m}}{\partial x_{j}}+\frac{\partial u^{(1)}_{i}}{\partial x_{m}}\frac{\partial u^{(2)}_{j}}{\partial x_{m}}+\frac{\partial u^{(2)}_{i}}{\partial x_{m}}\frac{\partial u^{(1)}_{j}}{\partial x_{m}}\\ &\quad\quad\quad\quad\quad\quad\quad\quad\quad+\frac{\partial u^{(1)}_{i}}{\partial x_{m}}\frac{\partial u^{(2)}_{m}}{\partial x_{j}}+\frac{\partial u^{(2)}_{i}}{\partial x_{m}}\frac{\partial u^{(1)}_{m}}{\partial x_{j}}\Bigg{)}\frac{\partial u^{(0)}_{i}}{\partial x_{j}}.\end{split}

Linearized traction-to-displacement map. Let $v$ be the solution to the initial boundary value problem with Neumann boundary value

(35)

\begin{split}&\rho\frac{\partial^{2}v}{\partial t^{2}}-\nabla\cdot S^{L}(x,v)=0,\quad(t,x)\in(0,T)\times\Omega,\\ &\nu\cdot S^{L}(x,v)=g,\quad(t,x)\in(0,T)\times\partial\Omega,\\ &v(0,x)=\frac{\partial}{\partial t}v(0,x)=0,\quad x\in\Omega,\end{split}

where $g$ is supported away from $\{t=0\}$ . Denote the Neumann-to-Dirichlet map

\Lambda^{{}^{\prime}}_{N}:g\mapsto v|_{(0,T)\times\partial\Omega}.

It is clear to see that

\Lambda^{{}^{\prime}}_{N}=(\Lambda^{{}^{\prime}}_{D})^{-1}.

Let $v^{(j)},j=1,2$ be solution to (35) with $g=g^{(1)}$ , and $v^{(12)}$ be solution to the following initial boundary value problem

(36)

\begin{split}&\rho\frac{\partial^{2}v^{(12)}}{\partial t^{2}}-\nabla\cdot S^{L}(x,u^{(12)})=\nabla\cdot G(v^{(1)},v^{(2)}),\quad(t,x)\in(0,T)\times\Omega,\\ &\nu\cdot S^{L}(x,v^{(12)})=-\nu\cdot G(v^{(1)},v^{(2)}),\quad(t,x)\in(0,T)\times\partial\Omega,\\ &v^{(12)}(0,x)=\frac{\partial}{\partial t}v^{(12)}(0,x)=0,\quad x\in\Omega.\end{split}

Denote then

\Lambda^{{}^{\prime\prime}}_{N}:(g^{(1)},g^{(2)})\mapsto v^{(12)}|_{(0,T)\times\partial\Omega}.

Assume $v^{(0)}$ solves the initial boundary value problem for the backward elastic wave equation

(37)

\begin{split}&\rho\frac{\partial^{2}}{\partial t^{2}}v^{(0)}-\nabla\cdot S^{L}(x,v^{(0)})=0,\quad(t,x)\in(0,T)\times\Omega,\\ &\nu\cdot S^{L}(x,v^{(0)})(t,x)=g^{(0)},\quad(t,x)\in(0,T)\times\partial\Omega,\\ &v^{(0)}(T,x)=\frac{\partial}{\partial t}v^{(0)}(T,x)=0,\quad x\in\Omega.\end{split}

Using integration by parts, we calculate

(38)

\begin{split}&\int_{0}^{T}\int_{\partial\Omega}\Lambda^{{}^{\prime\prime}}_{N}(g^{(1)},g^{(2)})\cdot g^{(0)}\mathrm{d}S\mathrm{d}t\\ =&\int_{0}^{T}\int_{\partial\Omega}v^{(12)}\cdot(\nu\cdot S^{L}(x,v^{(0)}))\mathrm{d}S\mathrm{d}t\\ =&\int_{0}^{T}\int_{\Omega}(\nabla\cdot S^{L}(x,v^{(0)}))\cdot v^{(12)}\mathrm{d}x\mathrm{d}t+\int_{0}^{T}\int_{\Omega}S^{L}(x,v^{(0)}):\nabla v^{(12)}\mathrm{d}x\mathrm{d}t\\ =&\int_{0}^{T}\int_{\Omega}\rho\frac{\partial^{2}v^{(0)}}{\partial t^{2}}\cdot v^{(12)}\mathrm{d}x\mathrm{d}t-\int_{0}^{T}\int_{\Omega}v^{(0)}\cdot(\nabla\cdot S^{L}(x,v^{(12)}))\mathrm{d}x\mathrm{d}t\\ &\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\int_{0}^{T}\int_{\partial\Omega}v^{(0)}\cdot(\nu\cdot S^{L}(x,v^{(12)})\mathrm{d}S\mathrm{d}t\\ =&\int_{0}^{T}\int_{\Omega}v^{(0)}\cdot\left(\rho\frac{\partial^{2}v^{(12)}}{\partial t^{2}}-\nabla\cdot S^{L}(x,v^{(12)})\right)\mathrm{d}x\mathrm{d}t-\int_{0}^{T}\int_{\partial\Omega}v^{(0)}\cdot(\nu\cdot G(v^{(1)},v^{(2)}))\mathrm{d}S\mathrm{d}t\\ =&\int_{0}^{T}\int_{\Omega}v^{(0)}\cdot(\nabla\cdot G(v^{(1)},v^{(2)}))\mathrm{d}x\mathrm{d}t-\int_{0}^{T}\int_{\partial\Omega}v^{(0)}\cdot(\nu\cdot G(v^{(1)},v^{(2)}))\mathrm{d}S\mathrm{d}t\\ =&-\int_{0}^{T}\int_{\Omega}\mathcal{G}(\nabla v^{(1)},\nabla v^{(2)},\nabla v^{(0)})\mathrm{d}x\mathrm{d}t,\end{split}

Note that

\begin{split}&\int_{0}^{T}\int_{\partial\Omega}\Lambda^{{}^{\prime\prime}}_{N}(g^{(1)},g^{(2)})\cdot g^{(0)}\mathrm{d}S\mathrm{d}t\\ =&-\int_{0}^{T}\int_{\Omega}\mathcal{G}(\nabla v^{(1)},\nabla v^{(2)},\nabla v^{(0)})\mathrm{d}x\mathrm{d}t\\ =&-\int_{0}^{T}\int_{\partial\Omega}\Lambda^{{}^{\prime\prime}}_{D}(v^{(1)},v^{(2)})\cdot v^{(0)}\mathrm{d}S\mathrm{d}t\\ =&-\int_{0}^{T}\int_{\partial\Omega}\Lambda^{{}^{\prime\prime}}_{D}(\Lambda^{{}^{\prime}}_{N}(g^{(1)}),\Lambda^{{}^{\prime}}_{N}(g^{(2)}))\cdot\Lambda^{{}^{\prime}}_{N}(g^{(0)})\mathrm{d}S\mathrm{d}t.\end{split}

Therefore, for the rest of the paper, we only need to consider the problem of recovering the parameters $\lambda,\mu,\rho,\mathscr{A},\mathscr{B},\mathscr{C}$ from $\Lambda^{{}^{\prime}}_{D}$ and $\Lambda^{{}^{\prime\prime}}_{N}$ .

Now assume that $\Lambda$ is the displacement-to-traction map associated with $\lambda,\mu,\rho,\mathscr{A},\mathscr{B},\mathscr{C}$ , and $\widetilde{\Lambda}$ is the displacement-to-traction map associated with $\widetilde{\lambda},\widetilde{\mu},\widetilde{\rho},\widetilde{\mathscr{A}},\widetilde{\mathscr{B}},\widetilde{\mathscr{C}}$ . Our goal is to prove that $\Lambda=\widetilde{\Lambda}$ implies

(\lambda,\mu,\rho,\mathscr{A},\mathscr{B},\mathscr{C})=(\widetilde{\lambda},\widetilde{\mu},\widetilde{\rho},\widetilde{\mathscr{A}},\widetilde{\mathscr{B}},\widetilde{\mathscr{C}}).

4.2. Uniqueness of the wavespeeds

Now assume $\Lambda^{\prime}_{D}=\widetilde{\Lambda}^{\prime}_{D}$ (or equivalently $\Lambda^{\mathrm{lin}}=\widetilde{\Lambda}^{\mathrm{lin}}$ ) and denote

g_{P}=c_{P}^{-2}\mathrm{d}s^{2}=\frac{\rho}{\lambda+2\mu}\mathrm{d}s^{2},\quad g_{S}=c_{S}^{-2}\mathrm{d}s^{2}=\frac{\rho}{\mu}\mathrm{d}s^{2},

and similarly

\widetilde{g}_{P}=\widetilde{c}_{P}^{-2}\mathrm{d}s^{2}=\frac{\widetilde{\rho}}{\widetilde{\lambda}+2\widetilde{\mu}}\mathrm{d}s^{2},\quad\widetilde{g}_{S}=\widetilde{c}_{S}^{-2}\mathrm{d}s^{2}=\frac{\widetilde{\rho}}{\widetilde{\mu}}\mathrm{d}s^{2}.

For the uniqueness of the two wavespeeds, we summarize the results in [29], [35] in the following proposition.

Proposition 1.

Assume that $(\Omega,g_{\bullet})$ and $(\Omega,\widetilde{g}_{\bullet})$ , where $\bullet=P,S$ , satisfy either of the following two conditions

(1)

$(\Omega,g_{\bullet}/\widetilde{g}_{\bullet})$ is simple;
(2)

$(\Omega,g_{\bullet}/\widetilde{g}_{\bullet})$ satisfies the strictly convex foliation condition.

Then $\Lambda^{\mathrm{lin}}=\widetilde{\Lambda}^{\mathrm{lin}}$ implies that

(39)

c_{P}=\widetilde{c}_{P},\quad c_{S}=\widetilde{c}_{S},\text{ or equivalently, }\quad\frac{\mu}{\rho}=\frac{\widetilde{\mu}}{\widetilde{\rho}},\quad\frac{\lambda}{\rho}=\frac{\widetilde{\lambda}}{\widetilde{\rho}},\quad\text{in }\overline{\Omega}.

From now on we assume the equality (39) to hold, and therefore $g_{\bullet}=\widetilde{g}_{\bullet}$ for $\bullet=S,P$ . We only need to work with the metrics $g_{\bullet}=c_{\bullet}^{-2}\mathrm{d}s^{2}$ in the following.

4.3. An integral identity

Fix $k$ large enough. By previous section, we can construct $u_{\varrho}$ and $\widetilde{u}_{\varrho}$ such that

(40)

\|\mathcal{L}_{\lambda,\mu,\rho}u_{\varrho}\|_{H^{k}}=\mathcal{O}(\varrho^{-K}),\quad\|\mathcal{L}_{\widetilde{\lambda},\widetilde{\mu},\widetilde{\rho}}\widetilde{u}_{\varrho}\|_{H^{k}}=\mathcal{O}(\varrho^{-K}),

where $u_{\varrho}=u_{\varrho}^{\mathrm{inc}}+u_{\varrho}^{\mathrm{ref}}$ , $\widetilde{u}_{\varrho}=\widetilde{u}_{\varrho}^{\mathrm{inc}}+\widetilde{u}_{\varrho}^{\mathrm{ref}}$ , and $u^{\mathrm{inc}}_{\varrho}$ and $\widetilde{u}^{\mathrm{inc}}_{\varrho}$ are Gaussian beam solutions concentrating near a null-geodesic $\vartheta$ in $(M,-\mathrm{d}t^{2}+c_{\bullet}^{-2}\mathrm{d}s^{2})$ . Furthermore we have

\|\mathcal{N}_{\lambda,\mu}u_{\varrho}\|_{H^{k}((0,T)\times\partial\Omega\setminus R_{0})}=\mathcal{O}(\varrho^{-K}),\quad\|\mathcal{N}_{\widetilde{\lambda},\widetilde{\mu}}\widetilde{u}_{\varrho}\|_{H^{k}((0,T)\times\partial\Omega\setminus R_{0})}=\mathcal{O}(\varrho^{-K}).

Since $\Lambda^{\mathrm{lin}}=\widetilde{\Lambda}^{\mathrm{lin}}$ , we know that jets of $(\lambda,\mu,\rho)$ and $(\widetilde{\lambda},\widetilde{\mu},\widetilde{\rho})$ at $\partial\Omega$ are equal (cf. [28]). Therefore we can denote

\mathcal{N}=\mathcal{N}_{\lambda,\mu}=\mathcal{N}_{\widetilde{\lambda},\widetilde{\mu}}.

Also we can extend $(\lambda,\mu,\rho)$ and $(\widetilde{\lambda},\widetilde{\mu},\widetilde{\rho})$ smoothly to $\widetilde{\Omega}$ such that $(\lambda,\mu,\rho)=(\widetilde{\lambda},\widetilde{\mu},\widetilde{\rho})$ on $\widetilde{\Omega}\setminus\Omega$ . By the constructions in Section 2, we can arrange the boundary values of $u_{\varrho}^{\mathrm{inc}}$ and $\widetilde{u}_{\varrho}^{\mathrm{inc}}$ such that

(41)

\|\mathcal{N}(u_{\varrho}-\widetilde{u}_{\varrho})\|_{H^{k}(R_{0})}=\|\mathcal{N}(u^{\mathrm{inc}}_{\varrho}-\widetilde{u}^{\mathrm{inc}}_{\varrho})\|_{H^{k}(R_{0})}\leq C\varrho^{-K}.

We refer to [16] for more details.

Let $f_{\varrho}=u_{\varrho}|_{(0,T)\times\partial\Omega}$ . Then one can construct solutions $u$ and $\widetilde{u}$ to

(42)

\begin{split}\mathcal{L}_{\lambda,\mu,\rho}u&=0,\quad\text{on }(0,T)\times\Omega,\\ \mathcal{N}u&=f_{\varrho},\quad\text{on }(0,T)\times\partial\Omega,\\ u(0,x)=\partial_{t}u(0,x)&=0,\quad\text{for }x\in\Omega.\end{split}

and

(43)

\begin{split}\mathcal{L}_{\widetilde{\lambda},\widetilde{\mu},\widetilde{\rho}}\widetilde{u}&=0,\quad\text{on }(0,T)\times\Omega,\\ \mathcal{N}\widetilde{u}&=f_{\varrho},\quad\text{on }(0,T)\times\partial\Omega,\\ \widetilde{u}(0,x)=\partial_{t}\widetilde{u}(0,x)&=0,\quad\text{for }x\in\Omega.\end{split}

To be more precise, we construct above solutions such that

u=u_{\varrho}+R_{\varrho},\quad\widetilde{u}=\widetilde{u}_{\varrho}+\widetilde{R}_{\varrho},

where $R_{\varrho}$ and $\widetilde{R}_{\varrho}$ satisfy

(44)

\begin{split}\mathcal{L}_{\lambda,\mu,\rho}R_{\varrho}&=-\mathcal{L}_{\lambda,\mu,\rho}u_{\varrho},\quad\text{on }(0,T)\times\Omega,\\ \mathcal{N}R&=0,\quad\text{on }(0,T)\times\partial\Omega,\\ R(0,x)=\partial_{t}R(0,x)&=0,\quad\text{for }x\in\Omega,\end{split}

and

(45)

\begin{split}\mathcal{L}_{\widetilde{\lambda},\widetilde{\mu},\widetilde{\rho}}\widetilde{R}_{\varrho}&=-\mathcal{L}_{\widetilde{\lambda},\widetilde{\mu},\widetilde{\rho}}\widetilde{u}_{\varrho},\quad\text{on }(0,T)\times\Omega,\\ \mathcal{N}\widetilde{R}&=(u_{\varrho}-\widetilde{u}_{\varrho})|_{(0,T)\times\partial\Omega},\quad\text{on }(0,T)\times\partial\Omega,\\ \widetilde{R}(0,x)=\partial_{t}\widetilde{R}(0,x)&=0,\quad\text{for }x\in\Omega.\end{split}

By (40) and (41), one can obtain

\|R_{\varrho}\|_{H^{k+1}},\|\widetilde{R}_{\varrho}\|_{H^{k+1}}\leq C\varrho^{-K}.

using standard theory for linear hyperbolic systems. Take $K$ large enough, then if $k>2$ we can use Sobolev imbedding to have

\|R_{\varrho}\|_{W^{1,3}},\|\widetilde{R}_{\varrho}\|_{W^{1,3}}=\mathcal{O}(\varrho^{-1/2}).

We note also that

\|u_{\varrho}\|_{W^{1,3}},\|\widetilde{u}_{\varrho}\|_{W^{1,3}}=\mathcal{O}(\varrho^{1/2}).

Let $\vartheta^{(j)}$ be a null-geodesic, and construct Gaussian beam solutions $u^{(j)}_{\kappa_{j}\varrho}$ and $\widetilde{u}^{(j)}_{\kappa_{j}\varrho}$ , whose incident waves $u^{(j),\mathrm{inc}}_{\kappa_{j}\varrho}$ and $\widetilde{u}^{(j),\mathrm{inc}}_{\kappa_{j}\varrho}$ concentrate near $\vartheta^{(j)}$ , where $\kappa_{j}$ is a constant. Then we construct $u^{(j)}$ and $\widetilde{u}^{(j)}$ , $j=1,2$ be solutions to (42) and (43) with $f_{\kappa_{j}\varrho}=u^{(j)}_{\kappa_{j}\varrho}|_{(0,T)\times\partial\Omega}$ . Also, let $\vartheta^{(0)}$ be a backward null-geodesic, and construct Gaussian beam solution $u^{(0)}_{\kappa_{0}\varrho}$ , $\widetilde{u}^{(0)}_{\kappa_{0}\varrho}$ , whose incident waves $u^{(0),\mathrm{inc}}_{\kappa_{0}\varrho}$ and $\widetilde{u}^{(0),\mathrm{inc}}_{\kappa_{0}\varrho}$ concentrate near $\vartheta^{(0)}$ . Then let $f^{(0)}_{\kappa_{0}\varrho}=u^{(0)}_{\kappa_{0}\varrho}|_{(0,T)\times\partial\Omega}$ and construct exact solutions $u^{(0)}$ , $\widetilde{u}^{(0)}$ to the backward elastic wave equations

(46)

\begin{split}\mathcal{L}_{\lambda,\mu,\rho}u^{(0)}&=0,\quad\text{on }(0,T)\times\Omega,\\ \mathcal{N}u^{(0)}&=f^{(0)}_{\kappa_{0}\varrho},\quad\text{on }(0,T)\times\partial\Omega,\\ u^{(0)}(T,x)=\partial_{t}u^{(0)}(T,x)&=0,\quad\text{for }x\in\Omega,\end{split}

and

(47)

\begin{split}\mathcal{L}_{\widetilde{\lambda},\widetilde{\mu},\widetilde{\rho}}\widetilde{u}^{(0)}&=0,\quad\text{on }(0,T)\times\Omega,\\ \mathcal{N}\widetilde{u}^{(0)}&=f^{(0)}_{\kappa_{0}\varrho},\quad\text{on }(0,T)\times\partial\Omega,\\ \widetilde{u}^{(0)}(T,x)=\partial_{t}\widetilde{u}^{(0)}(T,x)&=0,\quad\text{for }x\in\Omega.\end{split}

Then we carry out the second order linearization of traction-to-displacement map $\Lambda$ (and $\widetilde{\Lambda}$ ), the identity

\Lambda=\widetilde{\Lambda},

yields

\int_{0}^{T}\int_{\partial\Omega}\Lambda_{N}^{\prime\prime}(f^{(1)}_{\kappa_{1}\varrho},f^{(2)}_{\kappa_{2}\varrho})\cdot f^{(0)}_{\kappa_{0}\varrho}\mathrm{d}S\mathrm{d}t=\int_{0}^{T}\int_{\partial\Omega}\widetilde{\Lambda}_{N}^{\prime\prime}(f^{(1)}_{\kappa_{1}\varrho},f^{(2)}_{\kappa_{2}\varrho})\cdot f^{(0)}_{\kappa_{0}\varrho}\mathrm{d}S\mathrm{d}t,

and using (38) we obtain

(48)

\mathcal{I}:=\int_{0}^{T}\int_{\Omega}\mathcal{G}(\nabla u^{(1)},\nabla u^{(2)},\nabla u^{(0)})\,\mathrm{d}x\mathrm{d}t=\int_{0}^{T}\int_{\Omega}\widetilde{\mathcal{G}}(\nabla\widetilde{u}^{(1)},\nabla\widetilde{u}^{(2)},\nabla\widetilde{u}^{(0)})\,\mathrm{d}x\mathrm{d}t=:\widetilde{\mathcal{I}},

In the following, for the sake of simplicity we denote

\begin{split}u_{\kappa_{1}\varrho}^{(1),\mathrm{inc}}=e^{\mathrm{i}\kappa_{1}\varrho\varphi^{(1)}}\chi^{(1)}(\mathbf{a}^{(1)}+\mathcal{O}(\varrho^{-1})),\\ u_{\kappa_{2}\varrho}^{(2),\mathrm{inc}}=e^{\mathrm{i}\kappa_{2}\varrho\varphi^{(2)}}\chi^{(2)}(\mathbf{a}^{(2)}+\mathcal{O}(\varrho^{-1})),\\ u_{\kappa_{0}\varrho}^{(0),\mathrm{inc}}=e^{\mathrm{i}\kappa_{0}\varrho\varphi^{(0)}}\chi^{(0)}(\mathbf{a}^{(0)}+\mathcal{O}(\varrho^{-1})),\end{split}

where $\mathbf{a}^{(j)}:=\mathbf{a}_{0}^{(j)}$ . If the support of $u_{\kappa_{1}\varrho}^{(1)}u_{\kappa_{2}\varrho}^{(2)}u_{\kappa_{0}\varrho}^{(0)}$ is equal to the support of $u_{\kappa_{1}\varrho}^{(1),\mathrm{inc}}u_{\kappa_{2}\varrho}^{(2),\mathrm{inc}}u_{\kappa_{0}\varrho}^{(0),\mathrm{inc}}$ (which is the case in the following proof), we have (see [41] for more details)

(49)

\varrho^{-1}\frac{\mathrm{i}}{\kappa_{1}\kappa_{2}\kappa_{3}}\mathcal{I}=\rho^{2}\int_{0}^{T}\int_{\Omega}e^{\mathrm{i}\varrho S}\chi^{(1)}\chi^{(2)}\chi^{(0)}\mathcal{A}\mathrm{d}x\mathrm{d}t+\mathcal{O}(\varrho^{-1/2}),

where

(50)

S=\kappa_{1}\varphi^{(1)}+\kappa_{2}\varphi^{(2)}+\kappa_{0}\varphi^{(0)}

and

\mathcal{A}:=\mathcal{G}(\mathbf{a}^{(1)}\otimes\nabla\varphi^{(1)},\mathbf{a}^{(2)}\otimes\nabla\varphi^{(2)},\mathbf{a}^{(0)}\otimes\nabla\varphi^{(0)})

with

\begin{split}&\mathcal{G}(\mathbf{a}^{(1)}\otimes\nabla\varphi^{(1)},\mathbf{a}^{(2)}\otimes\nabla\varphi^{(2)},\mathbf{a}^{(0)}\otimes\nabla\varphi^{(0)})\\ =&\mathscr{B}[(\mathbf{a}^{(1)}\cdot\nabla\varphi^{(1)})(\mathbf{a}^{(2)}\cdot\nabla\varphi^{(0)})(\mathbf{a}^{(0)}\cdot\nabla\varphi^{(2)})+(\mathbf{a}^{(2)}\cdot\nabla\varphi^{(2)})(\mathbf{a}^{(1)}\cdot\nabla\varphi^{(0)})(\mathbf{a}^{(0)}\cdot\nabla\varphi^{(1)})\\ &\quad\quad\quad\quad+(\mathbf{a}^{(2)}\cdot\nabla\varphi^{(1)})(\mathbf{a}^{(1)}\cdot\nabla\varphi^{(2)})(\mathbf{a}^{(0)}\cdot\nabla\varphi^{(0)})]\\ &+\frac{\mathscr{A}}{4}\left((\mathbf{a}^{(2)}\cdot\nabla\varphi^{(1)})(\mathbf{a}^{(1)}\cdot\nabla\varphi^{(0)})(\mathbf{a}^{(0)}\cdot\nabla\varphi^{(2)})+(\mathbf{a}^{(1)}\cdot\nabla\varphi^{(2)})(\mathbf{a}^{(2)}\cdot\nabla\varphi^{(0)})(\mathbf{a}^{(0)}\cdot\nabla\varphi^{(1)})\right)\\ &+(\lambda+\mathscr{B})[(\mathbf{a}^{(1)}\cdot\nabla\varphi^{(1)})(\mathbf{a}^{(2)}\cdot\mathbf{a}^{(0)})(\nabla\varphi^{(2)}\cdot\nabla\varphi^{(0)})+(\mathbf{a}^{(2)}\cdot\nabla\varphi^{(2)})(\mathbf{a}^{(1)}\cdot\mathbf{a}^{(0)})(\nabla\varphi^{(1)}\cdot\nabla\varphi^{(0)})\\ &\quad\quad+(\mathbf{a}^{(1)}\cdot\mathbf{a}^{(2)})(\nabla\varphi^{(1)}\cdot\nabla\varphi^{(2)})(\mathbf{a}^{(0)}\cdot\nabla\varphi^{(0)})]+2\mathscr{C}(\mathbf{a}^{(1)}\cdot\nabla\varphi^{(1)})(\mathbf{a}^{(2)}\cdot\nabla\varphi^{(2)})(\mathbf{a}^{(0)}\cdot\nabla\varphi^{(0)})\\ &+(\mu+\frac{\mathscr{A}}{4})\Big{(}(\mathbf{a}^{(1)}\cdot\mathbf{a}^{(2)})(\nabla\varphi^{(1)}\cdot\mathbf{a}^{(0)})(\nabla\varphi^{(2)}\cdot\nabla\varphi^{(0)})+(\nabla\varphi^{(1)}\cdot\nabla\varphi^{(2)})(\mathbf{a}^{(1)}\cdot\mathbf{a}^{(0)})(\mathbf{a}^{(2)}\cdot\nabla\varphi^{(0)})\\ &\quad\quad+(\mathbf{a}^{(2)}\cdot\mathbf{a}^{(1)})(\nabla\varphi^{(2)}\cdot\mathbf{a}^{(0)})(\nabla\varphi^{(1)}\cdot\nabla\varphi^{(0)})+(\mathbf{a}^{(2)}\cdot\mathbf{a}^{(0)})(\mathbf{a}^{(1)}\cdot\nabla\varphi^{(0)})(\nabla\varphi^{(1)}\cdot\nabla\varphi^{(2)})\\ &\quad\quad+(\nabla\varphi^{(1)}\cdot\mathbf{a}^{(2)})(\nabla\varphi^{(2)}\cdot\nabla\varphi^{(0)})(\mathbf{a}^{(1)}\cdot\mathbf{a}^{(0)})+(\nabla\varphi^{(2)}\cdot\mathbf{a}^{(1)})(\nabla\varphi^{(1)}\cdot\nabla\varphi^{(0)})(\mathbf{a}^{(2)}\cdot\mathbf{a}^{(0)})\Big{)}.\end{split}

4.4. Construction of S-S-P waves

We first introduce the notation

L_{q}^{\bullet,\pm}M:=\{(\tau,\xi)\in T^{*}_{q}M,\,\tau^{2}=c_{\bullet}^{2}|\xi|^{2},\,\pm\tau>0\},

where $\bullet=P,S$ .

Fix a point $x_{0}\in\Omega$ and take $q=(\frac{T}{2},x_{0})\in M$ . Take $\vartheta^{(j)}:(t_{-}^{(j)},t_{+}^{(j)})\rightarrow\Omega$ be a null-geodesic in $((0,T)\times\Omega,\overline{g}_{S})$ for $j=1,2$ , and $\vartheta^{(0)}:(t_{+}^{(0)},t_{-}^{(0)})\rightarrow(0,T)\times\Omega$ be a backward null-geodesic in $((0,T)\times\Omega,\overline{g}_{P})$ , satisfying the following conditions:

(1)

$\vartheta^{(1)},\vartheta^{(2)},\vartheta^{(0)}$ intersect at $q$ , that is

$\vartheta^{(1)}\left(\frac{T}{2}\right)=\vartheta^{(2)}\left(\frac{T}{2}\right)=\vartheta^{(0)}\left(\frac{T}{2}\right)=p;$
(2)

denote $\zeta^{(j)}=\mathrm{d}\vartheta^{(j)}|_{q}$ , then $\zeta^{(1)},\zeta^{(2)}\in L^{S,+}_{q}M$ , $\zeta^{(0)}\in L^{P,-}_{q}M$ such that

$\kappa_{1}\zeta^{(1)}+\kappa_{2}\zeta^{(2)}+\kappa_{0}\zeta^{(0)}=0.$

with $\kappa_{1},\kappa_{2},\kappa_{0}$ not equal to zero at the same time, $\zeta^{(1)}$ and $\zeta^{(2)}$ are linearly independent.

Since $(\Omega,g_{\bullet})$ is non-trapping and $\partial\Omega$ is convex with respect to $g_{\bullet}$ , we can choose $\gamma^{(1)}$ and $\gamma^{(0)}$ such that

•

$\gamma^{(1)}([t_{-}^{(1)},\frac{T}{2}])$ , as a geodesic in $(\overline{\Omega},g_{S})$ , has no conjugate points on it;
•

$\gamma^{(0)}([\frac{T}{2},t_{-}^{(0)}])$ , as a geodesic in $(\overline{\Omega},g_{P})$ , has no conjugate points on it.

For $\zeta^{(1)},\zeta^{(0)}$ given, the vector $\zeta^{(2)}$ can be chosen in the following way (see also [11]). Take $\zeta^{(1)}\in L_{q}^{S,+}M$ , $\zeta^{(0)}\in L_{q}^{P,-}M$ . We write $\zeta^{(1)}=(\tau^{(1)},\xi^{(1)})$ , $\zeta^{(0)}=(\tau^{(0)},\xi^{(0)})$ , $\xi^{(1)},\xi^{(0)}\in\mathbb{R}^{3}$ , $|\xi^{(1)}|=|\xi^{(0)}|=1$ . Then we have

(\tau^{(1)})^{2}=c_{S}^{2}|\xi^{(1)}|^{2},\quad(\tau^{(0)})^{2}=c_{P}^{2}|\xi^{(0)}|^{2}.

Then, we have $\tau^{(1)}=c_{S}$ and $\tau^{(0)}=-c_{P}$ . Then we consider the vector $\zeta^{(2)}$ of the form $\zeta^{(2)}=a\zeta^{(1)}+b\zeta^{(0)}$ , $a,b\in\mathbb{R}$ . Without loss of generality, we take $b=1$ . In order for $\zeta^{(2)}$ to be in $L_{q}^{S,+}M$ and $\zeta^{(1)},\zeta^{(2)}$ are linearly independent, we need

(a\tau^{(1)}+\tau^{(0)})^{2}=c_{S}^{2}|a\xi^{(1)}+\xi^{(0)}|^{2}.

The above equation, by simple calculation, is equivalent to

2a\left(c_{S}^{2}\xi^{(1)}\cdot\xi^{(0)}+c_{S}c_{P}\right)=c_{P}^{2}-c_{S}^{2}.

Since $|c_{S}^{2}\xi^{(1)}\cdot\xi^{(0)}|\leq c_{S}^{2}<c_{S}c_{P}$ , the above equation always has a solution

a=\frac{c_{P}^{2}-c_{S}^{2}}{2(c_{S}^{2}\cos\psi+c_{S}c_{P})}\neq 0.

Then we get a vector $\zeta^{(2)}\in L_{p}^{S,+}M$ .

Assume $\vartheta$ is a null-geodesic with endpoints on $\partial M$ , we extend $\vartheta$ to a (collection of) broken null-geodesics $\Upsilon$ in the following inductive way. If $\vartheta([t_{0},t_{1}])\subset\Upsilon$ is a null-geodesic in $(\overline{M},\overline{g}_{\bullet})$ joining two boundary boundary points $\vartheta(t_{0})$ and $\vartheta(t_{1})$ with $t_{1}<T$ , then add the segments $\vartheta^{\mathrm{ref},P}$ and $\vartheta^{\mathrm{ref},S}$ to $\Upsilon$ , where

•

$\vartheta^{\mathrm{ref},P}:[t_{1},t_{2}^{P}]$ is a null-geodesic in $(\overline{M},\overline{g}_{P})$ connecting two boundary points $\vartheta^{\mathrm{ref},P}(t_{1})$ and $\vartheta^{\mathrm{ref},P}(t_{2}^{P})$ , with $\vartheta^{\mathrm{ref},P}(t_{1})=\vartheta(t_{1})$ and $\dot{\vartheta}^{\mathrm{ref},P}(t_{1})^{\flat}|_{T\partial M}=\dot{\vartheta}(t_{1})^{\flat}|_{T\partial M}$ ;
•

$\vartheta^{\mathrm{ref},S}:[t_{1},t_{2}^{S}]$ is a null-geodesic in $(\overline{M},\overline{g}_{S})$ connecting two boundary points $\vartheta^{\mathrm{ref},S}(t_{1})$ and $\vartheta^{\mathrm{ref},S}(t_{2}^{S})$ , with $\vartheta^{\mathrm{ref},S}(t_{1})=\vartheta(t_{1})$ and $\dot{\vartheta}^{\mathrm{ref},S}(t_{1})^{\flat}|_{T\partial M}=\dot{\vartheta}(t_{1})^{\flat}|_{T\partial M}$ .

Following the above procedure, we extend $\vartheta^{(1)}$ and $\vartheta^{(2)}$ to $\Upsilon^{(1)}$ and $\Upsilon^{(2)}$ , and extend $\vartheta^{(0)}$ backward to $\Upsilon^{(0)}$ . Notice that

\Upsilon^{(1)}([t_{-}^{(1)},\frac{T}{2}])=\vartheta^{(1)}([t_{-}^{(1)},\frac{T}{2}]),\quad\Upsilon^{(2)}([t_{-}^{(2)},\frac{T}{2}])=\vartheta^{(2)}([t_{-}^{(2)},\frac{T}{2}]),\quad\Upsilon^{(0)}([\frac{T}{2},t_{-}^{(0)}])=\vartheta^{(0)}([\frac{T}{2},t_{-}^{(0)}]).

Since $\gamma^{(1)}([t_{-}^{(1)},\frac{T}{2}])$ has no conjugate points, $\vartheta^{(1)}([t_{-}^{(1)},\frac{T}{2}])$ and $\vartheta^{(2)}([t_{-}^{(2)},\frac{T}{2}])$ intersect only at $p$ . Because $\gamma^{(0)}:[\frac{T}{2},t_{-}^{(0)}]\rightarrow\Omega$ is length-minimizing (w.r.t. the metric $g_{P}$ ) for any two points on it and $c_{P}>c_{S}$ , so $\Upsilon^{(j)}$ $(j=1,2)$ can not intersect $\vartheta^{(0)}((\frac{T}{2},t_{-}^{(0)}])$ . In conclusion, we know that $\Upsilon^{(1)},\Upsilon^{(2)},\Upsilon^{(0)}$ intersect only at the point $q$ , that is,

\Upsilon^{(1)}\cap\Upsilon^{(2)}\cap\Upsilon^{(0)}=\{q\}.

Now, construct $u^{(j)},\widetilde{u}^{(j)}$ , $j=1,2,0$ as in Section 4.3. Consider $u^{(j)}$ first, we can write

u^{(j)}=u^{(j)}_{\kappa_{j}\varrho}+R^{(j)}_{\kappa_{j}\varrho},\quad u^{(j)}_{\kappa_{j}\varrho}=u^{(j),\mathrm{inc}}_{\kappa_{j}\varrho}+u^{(j),\mathrm{ref}}_{\kappa_{j}\varrho},

for $j=1,2,0$ , where $u^{(j),\mathrm{inc}}_{\kappa_{j}\varrho}$ is a Gaussian beam solution concentrating near $\vartheta^{(j)}$ . We have similar expressions for $\widetilde{u}^{(j)}$ . We note that $u^{(1),\mathrm{inc}}_{\kappa_{j}\varrho},u^{(2),\mathrm{inc}}_{\kappa_{j}\varrho}$ represent two incident S-waves, $u^{(0),\mathrm{inc}}_{\kappa_{j}\varrho}$ represents an incident P-wave. By the above consideration, the solutions can be constructed such that the intersection of the supports of $u^{(1)}_{\kappa_{j}\varrho},u^{(2)}_{\kappa_{j}\varrho},u^{(0)}_{\kappa_{j}\varrho}$ is a small neighborhood of $p$ . By the constructions in Section 2, we have

\begin{split}\mathbf{a}^{(1)}|_{\vartheta^{(1)}}=&\det(Y_{S}^{(1)})^{-1/2}c_{S}^{-1/2}\rho^{-1/2}\mathbf{e}^{(1)},\\ \mathbf{a}^{(2)}|_{\vartheta^{(2)}}=&\det(Y_{S}^{(2)})^{-1/2}c_{S}^{-1/2}\rho^{-1/2}\mathbf{e}^{(2)},\\ \mathbf{a}^{(0)}|_{\vartheta^{(0)}}=&\det(Y_{P}^{(0)})^{-1/2}c_{P}^{-1/2}\rho^{-1/2}\nabla\varphi^{(0)},\ \end{split}

where $\mathbf{e}^{(j)}$ is a parallel vector field on $\gamma^{(j)}$ , which is normal to $\gamma^{(j)}$ , for $j=1,2$ . We can, without loss of generality, assume that

\nabla\varphi^{(j)}(q)=c_{S}^{-1}(q)\xi^{(j)},\quad\mathbf{e}^{(j)}(q)=c_{S}^{-1}(q)\alpha^{(j)},\quad j=1,2,\quad\nabla\varphi^{(0)}(q)=c_{P}^{-1}(q)\xi^{(0)},

with $\alpha^{(j)}\in\mathbb{R}^{3}$ , $|\alpha^{(j)}|=1$ , and $\alpha^{(j)}\perp\xi^{(j)}$ .

With the above choice of $\zeta^{(1)},\zeta^{(2)},\zeta^{(0)}$ we have (cf. [41, Lemma 4])

(1)

$S(q)=0$ ;
(2)

$(\partial_{t}S(q),\nabla S(q))=0$ ;
(3)

$\Im S(x)\geq cd(x,q)^{2}$ for $x$ in a neighborhood of $q$ , where $c>0$ is a constant,

where $S$ is defined in (50). Substituting these solutions into (48), and using the method of stationary phase, we end up with (cf. (49))

\varrho^{-1}\frac{\mathrm{i}}{\kappa_{1}\kappa_{2}\kappa_{3}}\mathcal{I}=c_{0}\mathcal{A}(q)+\mathcal{O}(\varrho^{-1/2})=c_{0}\widetilde{\mathcal{A}}(q)+\mathcal{O}(\varrho^{-1/2})=\varrho^{-1}\frac{\mathrm{i}}{\kappa_{1}\kappa_{2}\kappa_{3}}\widetilde{\mathcal{I}},

with some constant $c_{0}\neq 0$ . Letting $\varrho\rightarrow+\infty$ , we have

(51)

\mathcal{A}(q)=\widetilde{\mathcal{A}}(q)

We refer to [41] for more details.

4.5. Determination of $\rho^{-3/2}(\lambda+\mathscr{B})$ and $\rho^{-3/2}(4\mu+\mathscr{A})$

We assume $\alpha^{(1)}=\alpha^{(2)}=\alpha$ and $\alpha\perp\mathrm{span}\{\xi^{(2)},\xi^{(0)}\}$ . Since $\xi^{(1)}\in\mathrm{span}\{\xi^{(2)},\xi^{(0)}\}$ , $\alpha\cdot\xi^{(1)}=0$ . Now we compute

\begin{split}&\det(Y_{S}^{(1)})^{1/2}\det(Y_{S}^{(2)})^{1/2}\det(Y_{P}^{(0)})^{1/2}c_{P}^{5/2}c_{S}^{5}\rho^{3/2}\mathcal{A}(q)\\ =&\mathscr{B}(x_{0})[(\alpha\cdot\xi^{(1)})(\alpha\cdot\xi^{(0)})(\xi^{(0)}\cdot\xi^{(2)})+(\alpha\cdot\xi^{(2)})(\alpha\cdot\xi^{(0)})(\xi^{(0)}\cdot\xi^{(1)})\\ &\quad\quad\quad\quad+(\alpha\cdot\xi^{(1)})(\alpha\cdot\xi^{(2)})(\xi^{(0)}\cdot\xi^{(0)})]\\ &+\frac{\mathscr{A}}{4}(x_{0})\left((\alpha\cdot\xi^{(1)})(\alpha^{(1)}\cdot\xi^{(0)})(\xi^{(0)}\cdot\xi^{(2)})+(\alpha\cdot\xi^{(2)})(\alpha\cdot\xi^{(0)})(\xi^{(0)}\cdot\xi^{(1)})\right)\\ &+(\lambda+\mathscr{B})(x_{0})[(\alpha\cdot\xi^{(1)})(\alpha\cdot\xi^{(0)})(\xi^{(2)}\cdot\xi^{(0)})+(\alpha\cdot\xi^{(2)})(\alpha\cdot\xi^{(0)})(\xi^{(1)}\cdot\xi^{(0)})\\ &\quad\quad+(\alpha\cdot\alpha)(\xi^{(1)}\cdot\xi^{(2)})(\xi^{(0)}\cdot\xi^{(0)})]+2\mathscr{C}(x_{0})(\alpha\cdot\xi^{(1)})(\alpha\cdot\xi^{(2)})(\xi^{(0)}\cdot\xi^{(0)})\\ &+(\mu+\frac{\mathscr{A}}{4})(x_{0})\Big{(}(\alpha\cdot\alpha)(\xi^{(1)}\cdot\xi^{(0)})(\xi^{(2)}\cdot\xi^{(0)})+(\xi^{(1)}\cdot\xi^{(2)})(\alpha\cdot\xi^{(0)})(\alpha\cdot\xi^{(0)})\\ &\quad\quad+(\alpha\cdot\alpha)(\xi^{(2)}\cdot\xi^{(0)})(\xi^{(1)}\cdot\xi^{(0)})+(\alpha\cdot\xi^{(0)})(\alpha\cdot\xi^{(0)})(\xi^{(1)}\cdot\xi^{(2)})\\ &\quad\quad+(\xi^{(1)}\cdot\alpha)(\xi^{(2)}\cdot\xi^{(0)})(\alpha\cdot\xi^{(0)})+(\xi^{(2)}\cdot\alpha)(\xi^{(1)}\cdot\xi^{(0)})(\alpha\cdot\xi^{(0)})\Big{)}\\ =&(\lambda+\mathscr{B})(x_{0})\xi^{(1)}\cdot\xi^{(2)}+(2\mu+\frac{\mathscr{A}}{2})(x_{0})(\xi^{(1)}\cdot\xi^{(0)})(\xi^{(2)}\cdot\xi^{(0)}).\end{split}

Assume that

\xi^{(1)}\cdot\xi^{(0)}=\cos\psi,

then

\begin{split}&\xi^{(1)}\cdot\xi^{(2)}=\xi^{(1)}\cdot(a\xi^{(1)}+\xi^{(0)})=a|\xi^{(1)}|^{2}+\xi^{(1)}\cdot\xi^{(0)}=\frac{c_{P}^{2}-c_{S}^{2}}{2(c_{S}^{2}\cos\psi+c_{S}c_{P})}+\cos\psi,\quad\\ &\xi^{(2)}\cdot\xi^{(0)}=(a\xi^{(1)}+\xi^{(0)})\cdot\xi^{(0)}=a\xi^{(1)}\cdot\xi^{(0)}+|\xi^{(0)}|^{2}=1+\frac{c_{P}^{2}-c_{S}^{2}}{2(c_{S}^{2}\cos\psi+c_{S}c_{P})}\cos\psi.\end{split}

Thus we have

\begin{split}&\det(Y_{S}^{(1)})^{1/2}\det(Y_{S}^{(2)})^{1/2}\det(Y_{P}^{(0)})^{1/2}c_{P}^{5/2}c_{S}^{5}\rho^{3/2}\mathcal{A}(q)\\ =&(\lambda+\mathscr{B})\left[\frac{c_{P}^{2}-c_{S}^{2}}{2(c_{S}^{2}\cos\psi+c_{S}c_{P})}+\cos\psi\right]+(2\mu+\frac{\mathscr{A}}{2})\left[1+\frac{c_{P}^{2}-c_{S}^{2}}{2(c_{S}^{2}\cos\psi+c_{S}c_{P})}\cos\psi\right]\cos\psi.\end{split}

Similarly, if we replace $u^{(j)}$ by $\widetilde{u}^{(j)}$ in the above procedure, and notice that $u^{(j),\mathrm{inc}}$ and $\widetilde{u}^{(j),\mathrm{inc}}$ are Gaussian beam solutions along the same null-geodesic $\vartheta^{(j)}$ , we have

\begin{split}\widetilde{\mathbf{a}}^{(1)}|_{\vartheta^{(1)}}=&\det(Y_{S}^{(1)})^{-1/2}c_{S}^{-1/2}\widetilde{\rho}^{-1/2}\mathbf{e}^{(1)},\\ \widetilde{\mathbf{a}}^{(2)}|_{\vartheta^{(2)}}=&\det(Y_{S}^{(2)})^{-1/2}c_{S}^{-1/2}\widetilde{\rho}^{-1/2}\mathbf{e}^{(2)},\\ \widetilde{\mathbf{a}}^{(0)}|_{\vartheta^{(0)}}=&\det(Y_{P}^{(0)})^{-1/2}c_{P}^{-1/2}\widetilde{\rho}^{-1/2}\nabla\varphi^{(0)},\ \end{split}

Then have

\begin{split}&\det(Y_{S}^{(1)})^{1/2}\det(Y_{S}^{(2)})^{1/2}\det(Y_{P}^{(0)})^{1/2}c_{P}^{5/2}c_{S}^{5}\widetilde{\rho}^{3/2}\widetilde{\mathcal{A}}(q)\\ =&\left((\widetilde{\lambda}+\widetilde{\mathscr{B}})\left[\frac{c_{P}^{2}-c_{S}^{2}}{2(c_{S}^{2}\cos\psi+c_{S}c_{P})}+\cos\psi\right]+(2\widetilde{\mu}+\frac{\widetilde{\mathscr{A}}}{2})\left[1+\frac{c_{P}^{2}-c_{S}^{2}}{2(c_{S}^{2}\cos\psi+c_{S}c_{P})}\cos\psi\right]\cos\psi\right)(x_{0}).\end{split}

Consequently, we have the following identity from (51),

\begin{split}&\left((\lambda+\mathscr{B})\left[\frac{c_{P}^{2}-c_{S}^{2}}{2(c_{S}^{2}\cos\psi+c_{S}c_{P})}+\cos\psi\right]+(2\mu+\frac{\mathscr{A}}{2})\left[1+\frac{c_{P}^{2}-c_{S}^{2}}{2(c_{S}^{2}\cos\psi+c_{S}c_{P})}\cos\psi\right]\cos\psi\right)\rho^{-3/2}(x_{0})\\ =&\left((\widetilde{\lambda}+\widetilde{\mathscr{B}})\left[\frac{c_{P}^{2}-c_{S}^{2}}{2(c_{S}^{2}\cos\psi+c_{S}c_{P})}+\cos\psi\right]+(2\widetilde{\mu}+\frac{\widetilde{\mathscr{A}}}{2})\left[1+\frac{c_{P}^{2}-c_{S}^{2}}{2(c_{S}^{2}\cos\psi+c_{S}c_{P})}\cos\psi\right]\cos\psi\right)\widetilde{\rho}^{-3/2}(x_{0})\end{split}

By varying $\psi$ in the above identity, we end up with

(52)

\rho^{-3/2}(\lambda+\mathscr{B})(x_{0})=\widetilde{\rho}^{-3/2}(\widetilde{\lambda}+\widetilde{\mathscr{B}})(x_{0}),\quad\quad\rho^{-3/2}(4\mu+\mathscr{A})(x_{0})=\widetilde{\rho}^{-3/2}(4\widetilde{\mu}+\widetilde{\mathscr{A}})(x_{0}).

Since $x_{0}$ is an arbitrary point in $\Omega$ , and have

\rho^{-3/2}(4\mu+\mathscr{A})=\widetilde{\rho}^{-3/2}(4\widetilde{\mu}+\widetilde{\mathscr{A}}),\quad\text{in }\overline{\Omega}.

4.6. Determination of $\left(2\mu+2\mathscr{B}+\mathscr{A}\right)\rho^{-3/2}$

We can rescale $\xi^{(2)}$ such that $|\xi^{(2)}|=1$ . Assume $\alpha^{(1)},\alpha^{(2)}\in\mathrm{span}\{\xi^{(2)},\xi^{(0)}\}$ , $|\alpha^{(1)}|=|\alpha^{(2)}|=1$ and $\alpha^{(1)}\perp\xi^{(1)}$ , $\alpha^{(2)}\perp\xi^{(2)}$ . This means that $\xi^{(1)},\xi^{(2)},\xi^{(0)},\alpha^{(1)},\alpha^{(2)}$ are all in the same plane. Still denote

\xi^{(1)}\cdot\xi^{(0)}=\cos\psi,\quad\xi^{(1)}\cdot\xi^{(2)}=\cos\alpha.

Then we have

\begin{split}\alpha^{(2)}\cdot\xi^{(1)}&=\sin\alpha,\quad\alpha^{(1)}\cdot\xi^{(2)}=-\sin\alpha,\quad\alpha^{(1)}\cdot\xi^{(0)}=-\sin\psi,\\ &\alpha^{(2)}\cdot\xi^{(0)}=\sin(\alpha-\psi),\quad\alpha^{(1)}\cdot\alpha^{(2)}=\cos\alpha.\end{split}

Now we calculate

\begin{split}&\det(Y_{S}^{(1)})^{1/2}\det(Y_{S}^{(2)})^{1/2}\det(Y_{P}^{(0)})^{1/2}c_{P}^{5/2}c_{S}^{5}\rho^{3/2}\mathcal{A}(p)\\ =&\mathscr{B}(x_{0})(\alpha^{(2)}\cdot\xi^{(1)})(\alpha^{(1)}\cdot\xi^{(2)})(\xi^{(0)}\cdot\xi^{(0)})\\ &+\frac{\mathscr{A}}{4}(x_{0})\left((\alpha^{(2)}\cdot\xi^{(1)})(\alpha^{(1)}\cdot\xi^{(0)})(\xi^{(0)}\cdot\xi^{(2)})+(\alpha^{(1)}\cdot\xi^{(2)})(\alpha^{(2)}\cdot\xi^{(0)})(\xi^{(0)}\cdot\xi^{(1)})\right)\\ &+(\lambda+\mathscr{B})(x_{0})(\alpha^{(1)}\cdot\alpha^{(2)})(\xi^{(1)}\cdot\xi^{(2)})(\xi^{(0)}\cdot\xi^{(0)})\\ &+(\mu+\frac{\mathscr{A}}{4})(x_{0})\Big{(}(\alpha^{(1)}\cdot\alpha^{(2)})(\xi^{(1)}\cdot\xi^{(0)})(\xi^{(2)}\cdot\xi^{(0)})+(\xi^{(1)}\cdot\xi^{(2)})(\alpha^{(1)}\cdot\xi^{(0)})(\alpha^{(2)}\cdot\xi^{(0)})\\ &\quad\quad+(\alpha^{(1)}\cdot\alpha^{(2)})(\xi^{(2)}\cdot\xi^{(0)})(\xi^{(1)}\cdot\xi^{(0)})+(\alpha^{(2)}\cdot\xi^{(0)})(\alpha^{(1)}\cdot\xi^{(0)})(\xi^{(1)}\cdot\xi^{(2)})\\ &\quad\quad+(\xi^{(1)}\cdot\alpha^{(2)})(\xi^{(2)}\cdot\xi^{(0)})(\alpha^{(1)}\cdot\xi^{(0)})+(\xi^{(2)}\cdot\alpha^{(1)})(\xi^{(1)}\cdot\xi^{(0)})(\alpha^{(2)}\cdot\xi^{(0)})\Big{)}\\ =&-\mathscr{B}(x_{0})\sin^{2}\alpha+\frac{\mathscr{A}}{4}(x_{0})\left(-\sin\alpha\sin\psi\cos(\alpha-\psi)-\sin\alpha\sin(\alpha-\psi)\cos\psi\right)\\ &+(\lambda+\mathscr{B})(x_{0})\cos^{2}\alpha\\ &+(\mu+\frac{\mathscr{A}}{4})(x_{0})\Big{(}2\cos\alpha\cos\psi\cos(\alpha-\psi)-2\cos\alpha\sin\psi\sin(\alpha-\psi)\\ &\quad\quad\quad\quad\quad\quad\quad\quad-\sin\alpha\cos(\alpha-\psi)\sin\psi-\sin\alpha\cos\psi\sin(\alpha-\psi)\Big{)}\\ =&-\mathscr{B}(x_{0})\sin^{2}\alpha-\frac{\mathscr{A}}{4}(x_{0})\sin^{2}\alpha+(\lambda+\mathscr{B})(x_{0})\cos^{2}\alpha+(\mu+\frac{\mathscr{A}}{4})(x_{0})\left(2\cos^{2}\alpha-\sin^{2}\alpha\right)\\ =&\left(\lambda+2\mu+\mathscr{B}+\frac{\mathscr{A}}{2}\right)(x_{0})\cos^{2}\alpha-\left(\mu+\mathscr{B}+\frac{\mathscr{A}}{2}\right)(x_{0})\sin^{2}\alpha.\end{split}

Using similar argument as above, we obtain

(53)

\begin{split}&\left(\left(\lambda+2\mu+\mathscr{B}+\frac{\mathscr{A}}{2}\right)\cos^{2}\alpha-\left(\mu+\mathscr{B}+\frac{\mathscr{A}}{2}\right)\sin^{2}\alpha\right)\rho^{-3/2}(x_{0})\\ =&\left(\left(\widetilde{\lambda}+2\widetilde{\mu}+\widetilde{\mathscr{B}}+\frac{\widetilde{\mathscr{A}}}{2}\right)\cos^{2}\alpha-\left(\widetilde{\mu}+\widetilde{\mathscr{B}}+\frac{\widetilde{\mathscr{A}}}{2}\right)\sin^{2}\alpha\right)\widetilde{\rho}^{-3/2}(x_{0}).\end{split}

By (52), we already have

(54)

\left(\lambda+2\mu+\mathscr{B}+\frac{\mathscr{A}}{2}\right)\rho^{-3/2}(x_{0})=\left(\widetilde{\lambda}+2\widetilde{\mu}+\widetilde{\mathscr{B}}+\frac{\widetilde{\mathscr{A}}}{2}\right)\widetilde{\rho}^{-3/2}(x_{0}).

For $\alpha\neq 0$ , we obtain from (53) and (54)

(55)

\left(2\mu+2\mathscr{B}+\mathscr{A}\right)\rho^{-3/2}(x_{0})=\left(2\widetilde{\mu}+2\widetilde{\mathscr{B}}+\widetilde{\mathscr{A}}\right)\widetilde{\rho}^{-3/2}(x_{0}).

4.7. Determination of $\mathscr{C}$

Using (54) and (55), we can conclude that

\rho^{-3/2}(\lambda+\mu)=\widetilde{\rho}^{-3/2}(\widetilde{\lambda}+\widetilde{\mu})\quad\text{in }\overline{\Omega}.

Since $\frac{\mu}{\rho}=\frac{\widetilde{\mu}}{\widetilde{\rho}}$ and $\frac{\lambda}{\rho}=\frac{\widetilde{\lambda}}{\widetilde{\rho}}$ , we have

(\lambda,\mu,\rho)=(\widetilde{\lambda},\widetilde{\mu},\widetilde{\rho})\quad\text{in }\overline{\Omega}.

Using (52), we have

(\mathscr{A},\mathscr{B})=(\widetilde{\mathscr{A}},\widetilde{\mathscr{B}})\quad\text{in }\overline{\Omega}.

For the final step, we can apply the result in [41] to obtain

\mathscr{C}=\widetilde{\mathscr{C}}\quad\text{in }\overline{\Omega}

by inverting a Jacobi-weighted ray transform of the first kind on $(\Omega,g_{P})$ . See also [1] for more details.

Invertibility of this kind of weighted ray transform is proved in [14] under the no conjugate points assumption.

Now assume that $(\Omega,g_{P})$ satisfies the the foliation condition. Recall that $\partial\Omega$ is strictly convex with respect to $g_{P}$ . Let $\varpi\in C^{\infty}(\widetilde{\Omega})$ be a defining function of $\partial\Omega$ , such that $\varpi>0$ in $\Omega\setminus\partial\Omega$ , $\varphi<0$ in $\widetilde{\Omega}\setminus\overline{\Omega}$ , and $\varpi$ vanishes non-degenerately on $\partial\Omega$ . For any $p\in\partial\Omega$ , there exists a function $\tilde{x}\in C^{\infty}(\widetilde{\Omega})$ with $\tilde{x}(p)=0$ , $\mathrm{d}\tilde{x}(p)=-\mathrm{d}\varpi(p)$ , such that for $c>0$ small enough, $O_{p}=\{\tilde{x}>-c\}$ has no conjugate points. Then the above Jacobi-weighted ray transform is invertible on $O_{p}$ . A layer stripping scheme can be used to derive global invertibility [38, 27].

References

[1] S. Acosta, G. Uhlmann, and J. Zhai. Nonlinear ultrasound imaging modeled by a westervelt equation. SIAM Journal on Applied Mathematics, 82(2):408–426, 2022.
[2] T. Balehowsky, A. Kujanpää, M. Lassas, and T. Liimatainen. An inverse problem for the relativistic boltzmann equation. Communications in Mathematical Physics, 396(3):983–1049, 2022.
[3] G. Bao and H. Zhang. Sensitivity analysis of an inverse problem for the wave equation with caustics. Journal of the American Mathematical Society, 27(4):953–981, 2014.
[4] A. S. Barreto and P. Stefanov. Recovery of a general nonlinearity in the semilinear wave equation. arXiv preprint arXiv:2107.08513, 2021.
[5] M. Belishev and A. Katchalov. Boundary control and quasiphotons in the problem of reconstruction of a Riemannian manifold via dynamic data. Journal of Mathematical Sciences, 79(4):1172–1190, 1996.
[6] S. Bhattacharyya. Local uniqueness of the density from partial boundary data for isotropic elastodynamics. Inverse Problems, 34(12):125001, 2018.
[7] X. Chen, M. Lassas, L. Oksanen, and G. P. Paternain. Detection of hermitian connections in wave equations with cubic non-linearity. Journal of the European Mathematical Society, 24(7):2191–2232, 2021.
[8] X. Chen, M. Lassas, L. Oksanen, and G. P. Paternain. Inverse problem for the yang–mills equations. Communications in Mathematical Physics, 384:1187–1225, 2021.
[9] X. Chen, M. Lassas, L. Oksanen, and G. P. Paternain. Retrieving yang–mills–higgs fields in minkowski space from active local measurements. arXiv preprint arXiv:2204.12776, 2022.
[10] N. S. Dairbekov. Integral geometry problem for nontrapping manifolds. Inverse Problems, 22(2):431, 2006.
[11] M. de Hoop, G. Uhlmann, and Y. Wang. Nonlinear interaction of waves in elastodynamics and an inverse problem. Mathematische Annalen, 376(1-2):765–795, 2020.
[12] D. Dos Santos Ferreira, Y. Kurylev, M. Lassas, and M. Salo. The Calderón problem in transversally anisotropic geometries. Journal of the European Mathematical Society, 18(11):2579–2626, 2016.
[13] A. Feizmohammadi, J. Ilmavirta, Y. Kian, and L. Oksanen. Recovery of time dependent coefficients from boundary data for hyperbolic equations. J. Spectr. Theory, 11(3):1107–1143, 2021.
[14] A. Feizmohammadi and L. Oksanen. An inverse problem for a semi-linear elliptic equation in riemannian geometries. Journal of Differential Equations, 269(6):4683–4719, 2020.
[15] A. Feizmohammadi and L. Oksanen. Recovery of zeroth order coefficients in non-linear wave equations. Journal of the Institute of Mathematics of Jussieu, 21(2):367–393, 2022.
[16] P. Hintz, G. Uhlmann, and J. Zhai. The dirichlet-to-neumann map for a semilinear wave equation on lorentzian manifolds. Communications in Partial Differential Equations, 47(12):2363–2400, 2022.
[17] P. Hintz, G. Uhlmann, and J. Zhai. An inverse boundary value problem for a semilinear wave equation on lorentzian manifolds. International Mathematics Research Notices, 2022(17):13181–13211, 2022.
[18] A. Kachalov, Y. Kurylev, and M. Lassas. Inverse boundary spectral problems. CRC Press, 2001.
[19] A. Katchalov and Y. Kurylev. Multidimensional inverse problem with incomplete boundary spectral data. Communications in Partial Differential Equations, 23(1-2):27–59, 1998.
[20] Y. Kurylev, M. Lassas, L. Oksanen, and G. Uhlmann. Inverse problem for Einstein-scalar field equations. Duke Math. J., 171(16):3215–3282, 2022.
[21] Y. Kurylev, M. Lassas, and G. Uhlmann. Inverse problems for Lorentzian manifolds and non-linear hyperbolic equations. Inventiones mathematicae, 212(3):781–857, 2018.
[22] M. Lassas, T. Liimatainen, L. Potenciano-Machado, and T. Tyni. Uniqueness and stability of an inverse problem for a semi-linear wave equation. arXiv preprint arXiv:2006.13193, 2020.
[23] M. Lassas, T. Liimatainen, L. Potenciano-Machado, and T. Tyni. Stability estimates for inverse problems for semi-linear wave equations on lorentzian manifolds. arXiv preprint arXiv:2106.12257, 2021.
[24] M. Lassas, G. Uhlmann, and Y. Wang. Inverse problems for semilinear wave equations on Lorentzian manifolds. Communications in Mathematical Physics, 360(2):555–609, 2018.
[25] R. G. Muhometov and V. G. Romanov. On the problem of finding an isotropic riemannian metric in an n-dimensional space. In Doklady Akademii Nauk, volume 243, pages 41–44. Russian Academy of Sciences, 1978.
[26] G. P. Paternain, M. Salo, and G. Uhlmann. Invariant distributions, Beurling transforms and tensor tomography in higher dimensions. Mathematische Annalen, 363(1-2):305–362, 2015.
[27] G. P. Paternain, M. Salo, G. Uhlmann, and H. Zhou. The geodesic X-ray transform with matrix weights. American Journal of Mathematics, 141(6):1707–1750, 2019.
[28] L. Rachele. Boundary determination for an inverse problem in elastodynamics. Communications in Partial Differential Equations, 25(11-12):1951–1996, 2000.
[29] L. Rachele. An inverse problem in elastodynamics: uniqueness of the wave speeds in the interior. Journal of Differential Equations, 162(2):300–325, 2000.
[30] L. Rachele. Uniqueness of the density in an inverse problem for isotropic elastodynamics. Transactions of the American Mathematical Society, 355(12):4781–4806, 2003.
[31] A. Sá Barreto and P. Stefanov. Recovery of a cubic non-linearity in the wave equation in the weakly non-linear regime. Communications in Mathematical Physics, 392(1):25–53, 2022.
[32] A. Sá Barreto and Y. Wang. Singularities generated by the triple interaction of semilinear conormal waves. Analysis & PDE, 14(1):135–170, 2021.
[33] V. A. Sharafutdinov. Integral geometry of a tensor field on a manifold with upper-bounded curvature. Sibirskii Matematicheskii Zhurnal, 33(3):192–204, 1992.
[34] P. Stefanov, G. Uhlmann, and A. Vasy. Boundary rigidity with partial data. Journal of the American Mathematical Society, 29(2):299–332, 2016.
[35] P. Stefanov, G. Uhlmann, and A. Vasy. Local recovery of the compressional and shear speeds from the hyperbolic DN map. Inverse Problems, 34(1):014003, 2017.
[36] P. Stefanov, G. Uhlmann, and A. Vasy. Inverting the local geodesic x-ray transform on tensors. Journal d’Analyse Mathematique, 136(1):151–208, 2018.
[37] P. Stefanov, G. Uhlmann, and A. Vasy. The transmission problem in linear isotropic elasticity. Pure and Applied Analysis, 3(1):109–161, 2021.
[38] G. Uhlmann and A. Vasy. The inverse problem for the local geodesic ray transform. Inventiones mathematicae, 205(1):83–120, 2016.
[39] G. Uhlmann and Y. Wang. Determination of space-time structures from gravitational perturbations. Communications on Pure and Applied Mathematics, 73(6):1315–1367, 2020.
[40] G. Uhlmann and J. Zhai. Inverse problems for nonlinear hyperbolic equations. Discrete & Continuous Dynamical Systems: Series A, 41(1), 2021.
[41] G. Uhlmann and J. Zhai. On an inverse boundary value problem for a nonlinear elastic wave equation. Journal de Mathématiques Pures et Appliquées, 153:114–136, 2021.
[42] G. Uhlmann and Y. Zhang. Inverse boundary value problems for wave equations with quadratic nonlinearities. Journal of Differential Equations, 309:558–607, 2022.
[43] Y. Wang and T. Zhou. Inverse problems for quadratic derivative nonlinear wave equations. Communications in Partial Differential Equations, 44(11):1140–1158, 2019.

Determination of the density in a nonlinear elastic wave equation

Abstract.

Key words and phrases:

1. Introduction

Theorem 1.

2. Gaussian beams

2.1. Fermi coordinates

2.2. Gaussian beams

2.3. S-waves

2.4. P-waves

3. Gaussian beams with reflections

3.1. P- incident waves

Lemma 1.

Proof.

3.2. S- incident waves

3.3. Full Gaussian beam asymptotic solutions

4. Proof of the main theorem

4.1. Second order linearization of the displacement-to-traction map

4.2. Uniqueness of the wavespeeds

Proposition 1.

4.3. An integral identity

4.4. Construction of S-S-P waves

4.5. Determination of ρ−3/2​(λ+ℬ)\rho^{-3/2}(\lambda+\mathscr{B}) and ρ−3/2​(4​μ+𝒜)\rho^{-3/2}(4\mu+\mathscr{A})

4.6. Determination of (2​μ+2​ℬ+𝒜)​ρ−3/2\left(2\mu+2\mathscr{B}+\mathscr{A}\right)\rho^{-3/2}

4.7. Determination of 𝒞\mathscr{C}

References

4.5. Determination of $\rho^{-3/2}(\lambda+\mathscr{B})$ and $\rho^{-3/2}(4\mu+\mathscr{A})$

4.6. Determination of $\left(2\mu+2\mathscr{B}+\mathscr{A}\right)\rho^{-3/2}$

4.7. Determination of $\mathscr{C}$