Numerical approximation of the scattering amplitude in elasticity

We propose a numerical method to approximate the scattering amplitudes for the elasticity system with a non-constant matrix potential in dimensions $d=2$ and $3$. This requires to approximate first the scattering field, for some incident waves, which can be written as the solution of a suitable Lippmann-Schwinger equation. In this work we adapt the method introduced by G. Vainikko in \cite{V} to solve such equations when considering the Lam\'e operator. Convergence is proved for sufficiently smooth potentials. Implementation details and numerical examples are also given.


Introduction
The scattering of time-harmonic elastic waves in R d (d = 2, 3) is described through the solutions of the Lamé equation where u : R d → R d is the displacement vector, ω > 0 is the frequency and ∆ * u(x) = µ∆Iu(x) + (λ + µ)∇div u(x), with ∆I denoting the d × d diagonal matrix with the Laplace operator on the diagonal. The constants λ and µ are known as the Lamé parameters and depend on the underlying elastic properties. Throughout this paper we will assume that µ > 0 and 2µ + λ > 0 so that the operator ∆ * is strongly elliptic. The square matrix Q represents the action of some live loads on a bounded region. This can also represent inhomogeneities in the density of the elastic material ρ(x). In this case, Q takes the particular form Q = ω 2 (1 − ρ(x))I. Here we assume that each component of Q, q αβ (x), α, β = 1, ..., d, is real, compactly supported and belongs to L r (R d ) for some r > d/2.
When Q = 0 a solution of (1) is expressed as the sum of its compressional part u p and the shear part u s where with and They are solutions of the vectorial homogeneous Helmholtz equations ∆Iu p (x) + k 2 p u p (x) = 0, and ∆Iu s (x) + k 2 s u s (x) = 0, respectively. The values kp and ks are known as the associated speed of propagation of longitudinal waves (p-waves) and transverse waves (swaves) respectively. Note that compressional waves satisfy ∇ × u p = 0 while for the transverse ones we have div u s = 0. When Q = 0 solutions u of (1) can be interpreted as perturbations of the homogeneous ones, i.e. we write where ui is the incident wave, a solution of the homogeneous Lamé equation (i.e. with Q = 0), and v the scattered solution.
Unicity of scattered solutions v is obtained from radiation conditions when |x| → ∞. A natural choice is to assume that there are no reflections coming from infinity. As Q is compactly supported this condition coincides with the analogous imposed to solutions of ∆ * v(x) + ω 2 v(x) = 0 in an exterior domain. These are known as the outgoing Kupradze radiation conditions and they are equivalent to make both v p and v s (according to the decomposition (3)-(4)) satisfy the corresponding outgoing Sommerfeld radiation conditions for the Hemholtz equation, that is, (∂r − iks)v s = o(r −(d−1)/2 ), r = |x| → ∞.
It is known that, with certain conditions in Q, the system (1),(5) together with the outgoing Kupradze radiation conditions (6)-(7) has a unique solution u (see Theorem 2.1 and Proposition 3.1 of [4] and chapter 5 in [9]).
Problem (1), (6)- (7) and (5) is equivalent to the integral equation where Φ(x) is the fundamental tensor of the Lamé operator ∆ * + ω 2 I. This tensor is well-known for dimension d = 2, 3 (see [2] and [12]). We give its expression in section 7 below. The implicit integral equation for u in (8) is known as a Lippmann-Schwinger equation. A trigonometrical collocation method was proposed in [15] for the numerical approximation of this equation, when considering the scalar case and the fundamental solution of the Hemholtz equation, (we used this method in [6] to give approximations of the potential in the inverse quantum scattering problem). A similar collocation method, but without trigonometric polynomials, was proposed in [13]. In this paper we adapt the method in [15] to approximate the solutions of (8). Apart of the fact that we are dealing with a more complex vector problem, there are two main difficulties arising. The first one comes from the asymptotic bound for the Fourier coefficients of a suitable truncation of the Green tensor field associated to the Lamé operator (see Lemma 2 below). This is required to prove the convergence of the method. The second one comes from the fact that we slightly change the point of view in [15] where the numerical method gives an approximation of the product Qu, and therefore u can be computed only in the support of Q. As described in [15], the solution u can be extended to R d but the estimate of the numerical approximation blows up in the proximity of the boundary of the support of Q. In our approach we solve the problem directly on u and therefore we obtain a global approximation in R d .
As an application we approximate scattering amplitudes. Let us briefly describe what these objects are and their interest. As incident waves we usually consider plane waves either transverse (plane s-waves) with polarization vector ϕ ∈ S d−1 orthogonal to the wave direction θ ∈ S d−1 , or longitudinal plane waves (plane p-waves) Note that transversal waves depend on an extra angle parameter ϕ, not determined by θ, only for d = 3. From now on, we omit the dependence on ϕ when considering the specific case d = 2.
If up is the solution of (1) with up = u p i + vp, vp the scattered solution satisfying the outgoing Kupradze radiation condition, then vp = v p p + v s p and we have the asymptotic as |x| → ∞ (see section 2 of [4] or [5] where v p p,∞ and v s p,∞ are known as the longitudinal and transverse scattering amplitudes of up respectively. Note that the last argument in v p p,∞ (ω, θ, x/|x|) corresponds to a direction in S d−1 . This is interpreted as the direction in which the scattering amplitude is observed.
In a similar way, if us is the solution of (1) with us = u s i + vs, vs the scattered solution satisfying the outgoing Kupradze radiation condition then vs = v p s + v s s and we obtain the similar asymptotic to (11) and (12) and the longitudinal and transverse scattering amplitudes of us, v p s,∞ (ω, θ, ϕ, x/|x|), v s s,∞ (ω, θ, ϕ, x/|x|).
Scattering amplitudes can be easily measured in practice with seismographs situated far away from the support of the potential Q. Thus, a natural question is if we can derive information of the elastic material from these scattering amplitudes, [8], [10], [4], [5] and [3] for numerical results. This is known as the inverse scattering problem in elasticity, (analogous problems arise in quantum scattering, acoustic or electromagnetic obstacle scattering, [7]). In our case, the inverse problem would be to obtain information about the matrix Q.
To understand the problem let us mention how much information we can derive from these data. We can define 2 scattering amplitudes coming from longitudinal incident waves u p i . They depend on the variables (ω, θ, x/|x|) ∈ R + × S d−1 × S d−1 . This gives us 2d − 1 free parameters. As they are vector waves, each scattering amplitude provides d(2d − 1) parameters. Thus, we have 2d(2d−1) parameters that we can obtain from the longitudinal incident waves. Concerning transverse incident waves u s i we can define 2(d − 1) different scattering amplitudes. Each one depends now on the variables (ω, θ, ϕ, We have now 3d − 3 parameters for each vector wave, and therefore d(3d − 3) total parameters coming from each scattering amplitude. Thus, we have 2(d − 1)d(3d − 3) values that we can obtain from the transverse incident waves. Considering both longitudinal and transverse incident waves we sum up 2d(2d − 1) + 2(d − 1)d(3d − 3) scattering amplitudes that we can measure theoretically. In the particular case d = 2 we obtain 24 different scattering amplitudes, and 102 when d = 3. Now, if we want to determine Q from these scattering amplitudes we have to compute d 2 functions depending on d coordinates. In dimension d = 2 this requires information depending on 8 variables, while for d = 3, 27 variables are involved. We see that, in principle the inverse problem consisting in recovering the matrix Q from scattering data is overdetermined. Of course, some of these scattering amplitudes may not be independent and, in practice, we cannot measure all these parameters accurately but it is reasonable to think that we can combine part of this information to recover Q. Then, it is natural to restrict the scattering data. In the linear elasticity problem, it is most usual to deal with either fixed-angle or backscattering data. In the first case we fix the direction of the incident waves θ (which is no longer a free parameter) and in the latter, only scattering amplitudes with the last entry fixed as x/|x| = −θ are considered. Whether the scattering amplitudes allow to recover the potential Q in these cases is widely open. There are some results only for some particular situations as for instance if Q(x) = q(x)I, or when λ + µ = 0, where the Lamé operator can be reduced to a system of two independent vector Hemholtz equations, [4].
When reconstruction is not known, or difficult to obtain, it is sometimes possible to define approximations of Q from scattering amplitudes (see [5] for backscattering data and [4] for fixed angle data). Let us give an example in dimension d = 2. Note that scattering amplitudes for transversal incident waves will not depend on the angle ϕ in this case and we omit this variable in the following. We define the Fourier transform of the 2 × 2 matrix Q b (that we denote Q b ) from scattering amplitudes as follows. For ξ ∈ R 2 with ξ = −2ωθ, ω > 0 and θ ∈ S 1 , where vp,∞ (resp. vs,∞) is a linear combination of v p p,∞ (ωpp, θ, −θ) and v s p,∞ (ωps, θ, −θ), for certain energies ωpp and ωps (resp. v p s,∞ (ωpp, −θ) and v s s,∞ (ωps, −θ)), and {e1, e2} is the canonical basis in R 2 . The matrix Q b is known as the Born approximation of Q for backscattering data. There are not estimates on how good is the Born approximation in general. Most of the results involve asymptotic estimates for their Fourier tranforms, that allows to deduce that the Born approximation Q b and Q share the same singularities. In other words, we can use the Born approximation to recover the singularities of Q, [5].
It is worth mentioning that here we are restricting ourselves to the recovering of Q from scattering data. However, similar problems can be stated for more general situations as to recover the Lamé coefficients or other quantities depending on the elastic properties, [8]. This is a field where there are few results.
As we mentioned before, we apply the numerical approximation of the Lippmann-Scwinger equation (8) to simulate scattering amplitudes. In fact, as we show below (section 6) the scattering amplitudes can be recovered from the potential and the scattering solutions of (8) for the different incident waves. This allows to construct synthetic scattering amplitudes from numerical approximation of (8). In particular, we can test reconstruction algorithms and compute Born approximations that can be compared with the potential Q. This will be done in a forthcoming publication.
The rest of the paper is divided as follows: In section 2 we reduce the Lippmann-Schwinger equation (8) to a multiperiodic problem by localizing and periodic extension. In section 3 we introduce the finite trigonometric space used in the discretization. In section 4 we derive the trigonometric collocation method. Convergence is proved in section 5. In section 6 we show how to approximate the scattering amplitudes. In section 7 we estimate the Fourier coefficientes of the Green function involved in the multiperiodic version of (8). Finally, in section 8 we include some numerical examples that illustrate the convergence of the method and the approximation of some scattering amplitudes.

Reduction to a multiperiodic problem
We first write the Lippmann-Schwinger equation (8) where we have written the first term in the right hand side as a generic smooth function f : From now on we assume that, for some ρ > 0, the components of Q, Q αβ , and those of f = (f1, f2, ..., f d ) T satisfy, In particular this guarantees that the components of Q and f are continuous functions.
Given R > 2ρ, we define Problem (14) on GR is equivalent to the following: given f and the Solving (16) in x ∈ GR allows us to obtain v in GR . Once this is known, the solution in R d \GR can be obtained by simple integration using (8). This allows us to localize the problem to the bounded domain GR.
The idea now is to use (16) to approximate the solution v in the smaller ball x ∈ B(0, ρ) ⊂ GR. For those points, only the values of Φ in the ball B(0, 2ρ) are involved and therefore, changing Φ outside this ball does not affect to the solution v(x). This allows us to localize the Green tensor function to a compact support function in GR without changing the solution. We consider a smooth cutting of the form where ψ : [0, ∞) → R satisfy the conditions Analogously, we truncate f by considering f (x)ψ(|x|), that we still denote f to simplify.
A more drastic cutting using the characteristic function of the ball B(0, R) is also possible but this is more convenient as we comment below.
Once localized the problem we extend v, f , the components of the matrixes Q and K to R-periodic functions for which we still use same notation. Thus, we change (16) by a multiperiodic integral equation The smooth cutting of f and Φ(x) guarantees a smooth periodic extension from GR to R d . The solvability of (18) is obviously deduced from the one of (16). Moreover, once solved (18) we obtain the solution of (16) in B(0, ρ) and where f in (19) represents the original function, i.e. before cutting and periodizing.

Finite dimensional trigonometric space
The family of exponentials constitutes an orthonormal basis on L 2 (GR) with the norm We also introduce the space H η = H η (GR) which consists of dR−multiperiodic functions (distributions) having finite norm the Fourier coefficients of u.
We now introduce a finite dimensional approximation of H η . Let us consider h = 2R/N with N ∈ N and a mesh on GR with grid points jh, j ∈ Z d h and We also consider T h the finite dimensional subspace of trigonometric polynomials of the form Any v h ∈ T h can be represented either through the Fourier coefficients or the nodal values, where ϕ h,j (kh) = δ jk , more specifically For a given v h ∈ T h , the nodal valuesv h and the Fourier coefficientŝ v h are related by the discrete Fourier transform F h as follows, where, as usual, F h relates the sequence x(n) (n = (n1, n2, ..., n d )) with X(j) according to This definition coincides with the usual one in numerical codes (as MAT-LAB) up to a translation, since it considers instead j k = 0, ..., N − 1. This must be taken into account in the implementation. The orthogonal projection from H η to T h is defined by the formula (We need that v is continuous in order to be able to define

Trigonometric collocation method
The trigonometric collocation method to solve (18) reads, where f h = S h f and Note that here K and Q are d × d matrixes while w h is a column vector with each component w i h ∈ T h . The integral applies to each component of the integrand vector function.
From the numerical point of view there is a difficulty in the discrete formulation (21) since Qv h / ∈ (T h ) d and should be approximated. Therefore we modify the discrete formulation as follows, The operator K leaves invariant the subspace (T h ) d and (23) is consistent.
Since each one of the components in the matrix K is a periodic function the eigenfunctions of the associated convolution operator are known to be ϕj, defined in (20), while their eigenvalues K αβ (j) are the Fourier coefficients, i.e.
We can write (23) in matrix form too. Let .., f d h ) T as dN d -column vectors with the node values of the d components respectively. We also write the matrix form of the discrete Fourier transform as F h and observe that S h (Qv h ) can be interpreted as the vector obtained by the product of the components of the node values of Q. We obtain the matrix formulation, where Thus, the matrix form of the collocation method at the nodes reads, and at the Fourier coefficients,

Convergence of the trigonometric collocation method
In this section we prove a convergence result for the collocation method described above. The proof requires three preparatory lemmas.
The proof of this Lemma can be found in [15] or [14] (dimension two) and [11]. Lemma 2. The Fourier coefficients of the components of the matrix K defined in (17) satisfy the following estimate for some constant c > 0 that depends only on ω, λ, µ, R, ρ and d.
We prove this lemma in section 7 below.

Lemma 3.
Assume that all the components in Q satisfy Q αβ ∈ H η , with η > d/2 and let ε > 0 be such that η ≥ 2 − ε Then, where L((H η ) d , (H η ) d ) denotes the space of linear and continuous functions in (H η ) d (of course the constant cη,ε,R,Q also depends on the constant of the Lemma 2) .
Proof. As a consequence of Lemma 2 we have that K ∈ L((H η−(2−ε) ) d , (H η ) d ) for any η ∈ R and ε > 0. Assume that v ∈ (H η ) η . By Lemma 1, We have the following result Proof. This is a generalization of the analogous proof for the scalar case in [15]. From Lemmas 1 and 2 we have The existence and uniqueness for (18) is easily deduced from the inverse of the operator I − K(Q·) ∈ L((H η ) d , (H η ) d ). In fact this inverse exists since K(Q·) is compact and the homogeneous integral equation (18) with f = 0 has only the trivial solution, since if there were a non-zero solution v, it would be a solution of (8) with ui = 0 in B(0, ρ) and extend to R d .
Concerning the collocation equation (23), we write On the other hand by (29) we have then if we take h such that , then we can guarantee the existence of the inverse of the operator I − K(S h Q·) in L((H η ) d , (H η ) d ), (see Corollary 1.1.2 of [14]), and therefore the existence and uniqueness of solution v of (23) in (H η ) d . But if v is a solution of (23), then v is necessarily in (T h ) d .
We now obtain (30). First of all observe that from (18) and (23) we obtain , for sufficiently small h, it is enough to estimate the second hand term in this last expression. For the first two terms, by Lemma 1, we have, Now we estimate S h K − KS h as a linear operator in (H η ) d . Note that the projector operator P h commutes with K. To simplify we prove it only for the case d = 2. For w = j∈Z 2 (w 1 (j), w 2 (j)) T ϕj, we have In particular, (S h K − KS h ) = 0 on (T h ) d . Therefore, The two terms here are estimated in a similar way. For the first one, let ε be such that 0 < ε < 2, From the above inequality and (31) we obtain (30).
The trigonometric collocation provides a numerical method to approximate the solution of (18), which coincides with the solution of (8) in x ∈ B(0, ρ). Outside this ball we can approximate the solution of (8) by a suitable discretization of formula (19). For example, with the trapezoidal rule we obtain, The error can be easily estimated from Theorem 4, where c(x) = c|x| −(d−1)/2 , c independent of x and h and |x| −(d−1)/2 comes from the fundamental tensor of Lamé operator (see Section 7) and the asymptotic behaviour of H (1) ν (r) when r → ∞. Note that c(x) does not blow up as |x| → ρ, even if Φ(x − jh) becomes singular for some values of j. This is due to the fact that B(0, ρ) has compact support included in Q and therefore the terms Φ(x − jh)Q(jh) in (32) vanish for those j. This is in contrast with the method in [15], where the error bound is lost for |x| → ρ.

Numerical approximation of scattering data
As we said in the introduction we usually consider incident waves as plane waves either transverse (plane s-waves) or longitudinal (plane p-waves).
For an incident p-wave u p i in the form (10) we can define two different scattering amplitudes given in (11)- (12). These can be written as (see [4] or [5] where up is the solution of (8) with ui = u p i . Here Π x/|x| is the orthogonal projection onto the line defined by the vector x.
Analogously, for an incident s-wave in the form (9) we can define two different scattering amplitudes given in (13). These can be written as where now us is the solution of (8) with ui = u s i .

Fourier coefficients of the Green function
In this section we prove Lemma 2. We divide this section in two subsections where we consider separately the cases d = 2 and d = 3.

The case d = 2
The fundamental solution when d = 2 is given by (see [12]), where for x ∈ R 2 \{0} (identified with a matrix 1 × 2) the matrix J is given by x T x |x| 2 , and for each v > 0, the functions Φ1 and Φ2 are given by the Hankel function of first kind and order k.
On the other hand, we remind that Let a > 0 be a number to be determined later, j ∈ Z 2 and we assume that |j| −a ≤ min{1/2, ρ}. From now on C indicates a positive and universal constant which depends only on ω, Φn n = 1, 2, λ, µ, ρ, R, d and a.
We can write (39) We start by estimating the first integral in (39), essentially by using the size of the ball B(0, |j| −a ). From (37) and (38) Note that this first term satisfies the bound in (28) as soon as a ≥ 1.
To deal with the second integral in (39) we use that ϕj = − R 2 π 2 |j| 2 ∆ϕj and Green's formula, Note that we only have to prove that the three integrals in the right hand side of the second equality above can be bounded by C log |j|. The second one is obviously bounded since K αβ is C ∞ in 0 < |x| ≤ R. The same argument applies to the first integral in the annulus ρ ≤ |x| ≤ R, so that we can reduce ourselves to the subregion |j| −a < |x| < ρ. Thus, we only have to show that and where we have used that K αβ is equal to Φ αβ in B(0, ρ).

The case d = 3
The proof is analogous to the case d = 2 so that we omit some details. The fundamental solution in dimension d = 3 is given by (see [1]) and as in the case d = 2, It can be seen that where Ψ1 and Ψ2 are in C ∞ [0, ∞) and null when kp = ks.
From (46) and (47) we have Again C will indicate a universal constant, depending on the same parameters as in the case d = 2. Let j ∈ Z 3 , a > 0 to be chosen later and we assume that |j| −a ≤ min{1/2, ρ}.

Implementation and numerical experiments
The discretization of the Lippmann-Schwinger equation (14) reduces to the finite dimensional implicit system (25). The solution requires, 1. An approximation of the Fourier coefficientsK h .

2.
A solver for the linear system. Here we use the gmres routine in MATLAB.
To approximate the Fourier coefficientsK αβ we consider the trapezoidal rule. This can be interpreted as the discrete Fourier transform of K αβ in a suitable uniform mesh centered in the origin. Therefore we use the fft algorithm in MATLAB. We must be careful since the Green function is singular at the origin. This affects to the numerical approximation of the Fourier coefficients in (24), where the Green tensor appears in the integrand. In particular, if we use the trapezoidal rule to approximate the integral, we obtain a singular value for this integral. To avoid this problem we simply replace the value at the origin by zero. This does not affect significantly to the accuracy of the trapezoidal rule. For example, in dimension d = 2 the functions in the Green tensor have a logarithmic singularity at the origin. The error derived by our choice affects the quadrature formula only near x = 0. More precisely, it affects the ball centered at x = 0 and radius h (the mesh size) B(0, h). The real value of the integral in this ball is easily estimated by Thus, if we replace the integral by the value obtained with the trapezoidal rule and our choice, the error is basically the same as the one associated to the trapezoidal rule for smooth functions. Now we show some numerical experiments that illustrate the convergence of the numerical method in dimension d = 2 and provide numerical approximations of scattering data.

Experiment 1
In this first experiment we illustrate the convergence of the solution of the Lippmann-Schwinger equation by considering an analytical solution. We define the vector function Then, v = (∆ * + ω 2 I)g, f = v − g and Q = I solve the Lippmann-Schwinger equation (14). In table 1 we illustrate the convergence of the solution as N grows.
Experiment 2. Here we show an example of the scattering fields and amplitudes that can be obtained with the numerical method presented in this paper. Unfortunately we do not have analytical solutions in this case to compare the approximation.

Computational domain
Real approximation domain, Figure 1: Numerical results of experiment 2: Re (v s (x) · θ ⊥ ) corresponding to a transverse incident wave with incident angle θ = (cos π/4, sin π/4) and ω = 50, both in the computational domain (left) and the subdomain where the approximation to the continuous solution holds (right).
We fix the angle of an incident transverse wave θ = (cos π/4, sin π/4). In Figure 1  In the computational domain the solution is periodic, due to the trigonometric basis that we use. Of course, this is not the case for the inner subdomain where we approximate the real solution.
Now we compare both the real parts of the longitudinal Re (vs(x) · θ) and transverse Re( vs(x) · θ ⊥ ) solutions for transverse incident waves at different frequencies and the same angle θ = (cos π/4, sin π/4). Note that here the incident wave u s i = e iksθ·x θ ⊥ is complex and so it is the scattered field. The real part (resp. imaginary part) corresponds to the real part (resp. imaginary part) of the incident wave. In particular, we observe how larger incident frequencies provide more complicated solutions. A similar comparison for the longitudinal incident wave u s i = e iksθ·x θ is given in Figure 3.