Antoine Benoit
email: antoine.benoit@univ-littoral.fr

Weak well-posedness of hyperbolic boundary value problems in a strip: when instabilities do not meet the geometry

published or not. The documents may come

Weak well-posedness of hyperbolic boundary value problems in a strip: when instabilities do not meet the geometry

Introduction

The aim of this paper is to investigate the optimal loss of regularity for the following hyperbolic boundary value problem in the strip R d-1 × ]0, 1[:

         L(∂)u := ∂ t u + d j=1 A j ∂ j u = f (t, x , x d ) ∈ R t × R d-1 x × ]0, 1[ := Ω, B 0 u |x d =0 = g 0 (t, x ) ∈ R t × R d-1 x := ω, B 1 u |x d =1 = g 1 (t, x ) ∈ ω, u |t≤0 = 0 (x , x d ) ∈ R d-1 x × [0, 1] := Γ, (1)
where the coefficients A j ∈ M N ×N (R), B 0 ∈ M p×N (R) and B 1 ∈ M (N -p)×N (R) (we refer to Assumption 2.2 for the precise value of p).

By loss of regularity we mean that the (unique) solution u to (1) is not as regular (in terms of Sobolev spaces) as the source terms of (1), namely f , g 0 and g 1 . Such problems are referred as weakly well-posed problems, in contrast with problems for which the unique solution to (1) is as regular as the data of the system (referred as strongly well-posed problems).

Before to turn to a description of the results of this article let us recall what is known in the more classical geometry of the half-space.

The main result in the study of the strong well-posedness of the analogous of (1) in the half-space R d-1 × R + is due to Kreiss in [START_REF] Kreiss | Initial boundary value problems for hyperbolic systems[END_REF] where the author gives a full characterization of the boundary conditions leading to strong well-posedness (see also [START_REF] Sarason | On hyperbolic mixed problems[END_REF]). This characterization, the so-called uniform Kreiss-Lopatinskii condition, relies on the normal mode analysis of the boundary value problem in the half-space. Without enter into technical details (we refer to [START_REF] Benzoni-Gavage | Multidimensional hyperbolic partial differential equations[END_REF]- [START_REF] Chazarain | Introduction à la théorie des équations aux dérivées partielles linéaires[END_REF] for a complete study and to Section 2 for more informations) to the time/space half-space problem we associate its so-called resolvent frequency form (obtained after Laplace in time and Fourier in tangential space transforms):

d dx d u(x d ; ζ) = A (ζ) u(x d ; ζ) + A -1 d f (x d ; ζ) x d > 0, B 0 u |x d =0 (ζ) = g 0 (ζ), (2)
where ζ stands for the frequency parameter. The behaviour of the solution to the ordinary differential equation (2) strongly depends on the matrix A and [START_REF] Kreiss | Initial boundary value problems for hyperbolic systems[END_REF] is based on a very precise analysis of the possible structures of this matrix.

There are four different kinds of frequencies depending on the generalized eigenvalues of A (ζ), namely: elliptic, mixed, hyperbolic and glancing.

The problem (2) satisfies the uniform Kreiss-Lopatinskii condition and thus is strongly well-posed if this condition (see Assumption 2.3 for more details) holds for all possible kinds of frequency. However it has also been observed on many examples (see for example the ones of [START_REF] Benzoni-Gavage | Generic types and transitions in hyperbolic initial-boundary-value problems[END_REF]) that generically boundary conditions do not satisfy the uniform Kreiss-Lopatinskii condition.

A natural question is thus to determine what is the behaviour of the solution when the uniform Kreiss-Lopatinskii condition degenerates ?

The first possibility is that the problem generates a Hadamard instability and thus is ill-posed. However when the uniform Kreiss-Lopatinskii condition breaks down on one of the four kinds of frequency introduced above we may still have existence and uniqueness of the solution up to some losses of regularity.

Historically, this question has first been studied in [START_REF] Sablé-Tougeron | Existence pour un problème de l'élastodynamique Neumann non linéaire en dimension 2[END_REF] (in view of its applications in elastodynamics) in the case of a degeneracy of the uniform Kreiss-Lopatinskii condition for an elliptic frequency and a loss of one derivative on the boundary is shown.

Then the same question has been addressed in [START_REF] Coulombel | Stabilité multidimensionnelle d'interfaces dynamiques. Application aux transitions de phase liquide-vapeur[END_REF] and [START_REF] Coulombel | Well-posedness of hyperbolic initial boundary value problems[END_REF] for a degeneracy on mixed and on hyperbolic frequencies (in this context the main physical applications are phase transition and shock wave stability). For a mixed (resp. hyperbolic) degeneracy a loss of a half (resp. one) derivative in the interior appears in addition of the loss of one derivative on the boundary.

Once these estimations for the loss of regularity of the solution are stated an interesting remaining question is the one of their sharpness. This question has been widely studied in the litterature and we refer to [START_REF] Marcou | Rigorous weakly nonlinear geometric optics for surface waves[END_REF] for the elliptic case to [START_REF] Benoit | Geometric optics expansions for linear hyperbolic boundary value problems and optimality of energy estimates for surface waves[END_REF] for the mixed case and to [START_REF] Coulombel | Geometric optics expansions with amplification for hyperbolic boundary value problems: linear problems[END_REF] for the hyperbolic case (see also the formal analysis of [START_REF] Majda | A theory for spontaneous Mach stem formation in reacting shock fronts. I. The basic perturbation analysis[END_REF] and [START_REF] Majda | Nonlinear geometric optics for hyperbolic mixed problems[END_REF]) and the sharpness of all the estimates described so far is now established.

So that at present time, the characterization of weakly well-posed boundary value problems in the halfspace is achieved (up to the somewhat specific case of a degeneracy on a glancing frequency).

The common point of all these proofs is that they use WKB expansions (that is an approximate solution to (1) in the asymptotic of highly oscillating source terms) to saturate the known energy estimates in order to show their sharpness.

Concerning the strip geometry wa can show that the uniform Kreiss-Lopatinskii condition on each side of the strip is sufficient to ensure strong well-posedness (up to pay some exponential growth of the solution with respect to time (see [Bena] in which this growth is studied)). However, in the author's knowledge, the classification of the possible weakly well-posed problems in the strip has not been investigated at present time.

In this article we restrict our attention to a degeneracy of the uniform Kreiss-Lopatinskii condition on an elliptic frequency, while in [Benc] the degeneracy takes place on a hyperbolic frequency.

The main results of this article, namely Theorems 2.1 and 2.2, are the following:

• Theorem 2.1 establishes the weak well-posedness of (1) when one of the boundary conditions is responsible of a failure of the uniform Kreiss-Lopatinskii on an elliptic frequency1 .

• Theorem 2.2 shows that the energy estimate of Theorem 2.1, with a loss of one derivative on the boundary where the uniform Kreiss-Lopatinskii condition fails but any loss in the interior, is sharp.

Consequently in the elliptic degeneracy framework the situation in the strip geometry is not worse than the one in the half-space geometry. This result is not surprising and it is linked to the fact that, contrary to hyperbolic modes, elliptic modes can not propagate information/singularities for one side of the strip to the other. However, we believe that this result is interesting for several reasons.

Firstly because it gives a class of hyperbolic boundary value problems in the strip for which the loss of derivatives is at most one2 and in this framework one can expect to use iterative methods in order to treat non-linear problems.

Secondly, and in the author's opinion it is the most interesting point, because this non-amplification of singularities has to be compared with the result obtained in [Benc]. Indeed for a degeneracy of the uniform Kreiss-Lopatinskii condition on a hyperbolic frequency, [Benc] shows that the loss of derivatives of the solution has to be increasing with the final time of resolution of the system.

So the results of this paper combined with the one of [Benc] show that the behaviour of the solution can be totally different from one kind of degeneracy to the other (this is not really the case for weakly well-posed boundary value problems in the half-space for which the number of losses remain finite, no matter the time scale is). This is an example of the wide variety of possible phenomena for hyperbolic boundary value problems with several boundaries. This paper is organized as follows: in Section 2 we describe the main assumptions on (1) and state the main results, namely Theorem 2.1 and Theorem 2.2. Then Section 3 gives the main steps in the proof of Theorem 2.1 while Section 4 exposes the ones of Theorem 2.2.

Notations and statements

Some general notations

In the core of the text we use •

• as a short hand notation for: "there exists C > 0 such that • ≤ C•". When C depends on some (relevant) parameter λ we specify it by using the notation λ . For a, b ∈ R we write a, b := [a, b] ∩ Z.

For T > 0, we introduce:

Ω T := ]-∞, T ] × R d-1 x × ]0, 1[ and ω T := ]-∞, T ] × R d-1 x .
The frequency space used in this article is given by:

Ξ := ζ := (γ + iτ, η) ∈ C × R d-1 \ γ ≥ 0 \ {(0, 0)} ,
and for convenience we use the notation

Ξ 0 := Ξ ∩ {γ = 0}. For s ∈ R, γ > 0 and X ⊂ R t × R d-1
x × ]0, 1[ we define the following weighted (in time) Sobolev spaces by

H s γ (X) := u ∈ D (X) \ e -γt u ∈ H s (X) , and we set L 2 γ (X) := H 0 γ (X).

Main notations and assumptions

The first assumption on the strip problem (1) is an hyperbolicity assumption for L(∂)

Assumption 2.1 [Constant hyperbolicity]There exist M ∈ 1, N and λ 1 , ..., λ M analytic real valued functions on R d \ {0} and finally µ 1 , ..., µ M ≥ 1 such that

∀ξ ∈ S d-1 , det τ + d j=1 ξ j A j = M k=1 (τ + λ k (ξ)) µ k , with λ 1 < • • • < λ M and the eigenvalues λ k (ξ) are semi-simple.
Then a simplifying assumption is to restrict the study to non-characteristic boundary conditions.

Assumption 2.2 [Non-characteristic boundary]The matrix A d is invertible. Moreover p is the number of strictly positive eigenvalues of A d .

With Assumptions (2.1) and 2.2 in hand we define the resolvent matrix of (1) by:

∀ζ ∈ Ξ, A (ζ) := -A -1 d (γ + iτ )I + i d-1 j=1 η j A j . (3)
The result of [START_REF] Hersh | Mixed problems in several variables[END_REF] implies that for all ζ ∈ Ξ \ Ξ 0 , the stable (resp. unstable) subspace E s (ζ) (resp. E u (ζ)), that is the space generated by the eigenmodes associated to the generalized eigenvalues with strictly negative (resp. positive) real part, is well-defined and of constant dimension equals to p (resp. N -p).

Moreover the results of [START_REF] Kreiss | Initial boundary value problems for hyperbolic systems[END_REF] and [START_REF] Métivier | The block structure condition for symmetric hyperbolic systems[END_REF] show that the stable subspace E s (ζ) can be extended by continuity up to Ξ 0 . Of course in this extension some of the real parts may vanish. If such a situation occurs the associated frequency are called: mixed, hyperbolic or glancing (essentially depending on the number of eigenvalues that degenerate, we refer to [START_REF] Métivier | The block structure condition for symmetric hyperbolic systems[END_REF] for more details).

In this paper we are, however, interested in the case where all the eigenvalues keep a signed real part3 in this extension process, these frequencies are referred as elliptic frequencies. With more details we give the following definition: Definition 2.1 (Elliptic area) Under Assumptions 2.1 and 2.2 the elliptic area E is the set of ζ ∈ Ξ 0 such that there exist a neighbourhood V of ζ in Ξ and a matrix T ∈ GL N (C) regular on V such that:

∀ζ ∈ V , A (ζ) = T -1 (ζ)diag(A s (ζ), A u (ζ))T (ζ),
where A s , A u ∈ M p×p (C) are regular on V and satisfy that for all ζ ∈ V , Re(sp

(A s (ζ)) ⊂ ]-∞, 0[ and Re(sp(A u (ζ))) ⊂ ]0, +∞[.
In particular, for ζ ∈ E the decomposition

C N = E s (ζ) ⊕ E u (ζ), (4)
of [START_REF] Hersh | Mixed problems in several variables[END_REF] still holds and we denote by

Π s (ζ) (resp. Π u (ζ)) the projection on E s (ζ) (resp. E u (ζ))
with respect to (4).

Our last main assumption deals with the boundary conditions B 0 and B 1 :

Assumption 2.3

Main results

The main results of the article are the following5 :

Theorem 2.1 Under Assumptions 2.1,2.2 and 2.3 there exists γ > 0 such that for all γ ≥ γ and for all

f ∈ L 2 γ (Ω), g 0 ∈ H 1/2 γ (ω), g 1 ∈ L 2 γ (ω), the strip problem (1) admits a unique solution u ∈ L 2 γ (Ω) with traces u |x d =0 ∈ H -1/2 γ (ω) and u |x d =1 ∈ L 2 γ (ω)
. Moreover u satisfies the energy estimate: for all γ ≥ γ we have

γ u 2 L 2 γ (Ω) + γ u |x d =0 2 H -1/2 γ (ω) + u |x d =1 2 L 2 γ (ω) 1 γ f 2 L 2 γ (Ω) + 1 γ g 0 2 H 1/2 γ (ω) + g 1 2 L 2 γ (ω) . (5)

Remark

• If the side {x d = 1} also violates the uniform Kreiss-Lopatinskii condition in the elliptic area (see 2. of Assumption 2.3) then Theorem 2.1 still holds by assuming that g 1 ∈ H 1/2 γ (ω) and with the energy estimate: for all γ ≥ γ,

u 2 L 2 γ (Ω) + u |x d =0 2 H -1/2 γ (ω) + u |x d =1 2 H -1/2 γ (ω) 1 γ 2 f 2 L 2 γ (Ω) + g 0 2 H 1/2 γ (ω) + g 1 2 H 1/2 γ (ω) . (6)
• If the boundary condition(s) on which the Kreiss-Lopatinskii condition degenerates in the elliptic area is (are) homogeneous then a direct consequence of (5) (( 6)) is that we have the energy estimate without loss of derivative ∀γ ≥ γ, u 2

L 2 (Ω) 1 γ 2 f 2 L 2 (Ω) .
In such a situation we then have examples of hyperbolic boundary value problems in the strip without loss of derivatives even if the uniform Kreiss-Lopatinskii condition is not satisfied.

The following theorem shows that (5) is sharp. More precisely Theorem 2.2 Under Assumptions 2.1, 2.2 and 2.3, let s > 0, we assume that for all f ∈ L 2 (Ω T ), g 0 ∈ H s (ω T ), g 1 ∈ L 2 (ω T ) the strip problem (1) admits a solution u such that

u 2 L 2 (Ω T ) T f 2 L 2 (Ω T ) + g 0 2 H s (ω T ) + g 1 2 L 2 (ω T ) , (7)
then necessarily s ≥ 1/2.

3 Weak well-posedness of (1), proof of Theorem 2.1

The proof of Theorem 2.1 follows the classical proof to show the strong well-posedness of hyperbolic boundary value problems in the half-space already used in [START_REF] Kreiss | Initial boundary value problems for hyperbolic systems[END_REF]. It is a three steps proof:

1. we first show an a priori energy estimate for (1).

2. Then we show the existence of a weak solution.

3. At last we conclude to the uniqueness of the constructed weak solution, by showing that it satisfies the a priori energy estimate, thanks to a "weak=strong" lemma.

The fact that the same sketch of proof (even if, of course it involves many technical difficulties) also operates for weakly well-posed problems (see [START_REF] Sablé-Tougeron | Existence pour un problème de l'élastodynamique Neumann non linéaire en dimension 2[END_REF]- [START_REF] Coulombel | Stabilité multidimensionnelle d'interfaces dynamiques. Application aux transitions de phase liquide-vapeur[END_REF] and [START_REF] Coulombel | Well-posedness of hyperbolic initial boundary value problems[END_REF]) is quite interesting in its own and shows the robustness of the whole method.

In Paragraph 3.1 we give a simple proof of the a priori energy estimate (that is (5)) based on localization methods. The proof exposed here is (likely to be) the simplest possible one because we are not interested in the minimal exponential growth in time rate in the estimates.

Then Paragraph 3.2 contains some elements of proof concerning the existence and the "weak=strong" lemma.

A priori energy estimate for (1)

Let u be a smooth solution to (1). We introduce two cut-off functions χ 0 , χ

1 ∈ D(]-2, 2[ , [0, 1]) that satisfy χ 0 (x d ) = 1 if x d ∈ 0, 1 3 , 0 if x d ≥ 2 3 ,
and

χ 1 (x d ) = 0 if x d ≤ 1 3 , 1 if x d ∈ 2 3 , 1 .
Define χ 1/2 := 1 -χ 0 -χ 1 and for all j ∈ 0, 2 , u j/2 := χ j/2 u.

By construction u 1/2 solves the Cauchy problem for the hyperbolic operator L(∂) so that we have the standard energy estimate: there exists γ 1/2 > 0 such that for all γ ≥ γ 1/2

γ u 1/2 2 L 2 γ (Ω) 1 γ L(∂)u 1/2 2 L 2 γ (Ω) . (8)
Then u 0 (resp. u 1 ) solves the boundary value problem in the half-space

R d-1 × ]0, ∞[ (resp. R d-1 × ]-∞, 1[) so that from [ST88] (resp. [Kre70]
) we have the energy estimate: there exists γ 0 > 0 (resp. γ 1 > 0) such that for all γ ≥ γ 0 (resp. γ ≥ γ 1 ):

u 0 2 L 2 γ (Ω) + u 0 |x d =0 2 H -1/2 γ (ω) 1 γ 2 L(∂)u 0 2 L 2 γ (Ω) + B 0 u 0 |x d =0 2 H 1/2 γ (ω) , ( 9
) resp. γ u 1 2 L 2 γ (Ω) + u 1 |x d =1 2 L 2 γ (ω) 1 γ L(∂)u 1 2 L 2 γ (Ω) + B 1 u 1 |x d =1 2 L 2 γ (ω) . ( 10
)
However for all j ∈ 0, 2 we have L(∂)u j/2 L 2 γ (Ω)

A d f L 2 γ (Ω) + u L 2 γ (Ω)
. So that from the triangle inequality, the trace equalities u |x d =j = u j |x d =j (for j ∈ {0, 1}) combined with (8), ( 9) and (10) it turns out that for all γ ≥ max(γ 0 , γ 1/2 , γ 1 ):

γ u 2 L 2 γ (Ω) + γ u |x d =0 2 H -1/2 γ (ω) + u |x d =1 2 L 2 γ (ω) 1 γ f 2 L 2 γ (Ω) + 1 γ g 0 H 1/2 γ (ω) + g 1 2 L 2 γ (ω) + 1 γ u 2 L 2 γ (Ω) ,
so that up to choose γ large enough we can absorb the term 1 γ u 2 L 2 γ (Ω) from the right hand side into the left hand side and (5) follows for all regular solution to (1).

Existence and uniqueness

Once the a priori energy estimate for (1) is established, the existence and the uniqueness of the solution to (1) is proved exactly as exposed in [START_REF] Sablé-Tougeron | Existence pour un problème de l'élastodynamique Neumann non linéaire en dimension 2[END_REF] (in the simplified framework of constant coefficients in (1)).

A way to show the existence of a solution is to introduce a so-called dual problem of (1). More precisely we define:

           L * (∂)v := -∂ t + d j=1 A T j ∂ j v = f in Ω, C 0 v |x d =0 = g 0 for (t, x ) ∈ ω, C 1 v |x d =1 = g 1 for (t, x ) ∈ ω, v |t≤0 = 0 on Γ, (11)
where the boundary matrices are chosen in such a way that the following decompositions hold:

for ∈ {0, 1} , A d = C T N + M T B ,
where M , N , C ∈ M p×N (R) in order to have the Green formula (for all u, v regular enough):

L(∂)u, v -u, L * (∂)v = 1 =0 (-1) +1 N u |x d = , C v |x d = + B u |x d = , M v |x d = , (12)
where The end of the proof consists in showing that the dual problem (11) satisfies Assumptions 2.1, 2.2 and 2.3 in order to show that it satisfies the a priori energy estimate so that, essentially the Green formula (12) combined with Riesz representation theorem permits to show the existence of a weak solution to (1) 6 .

The fact that the dual problem (11) satisfies Assumptions 2.1 and 2.2 is readable from the definition of L * (∂). The fact that (11) satisfies Assumption 2.3 is a consequence of the definition of the Lopatinskii determinant for the dual problem (we refer to Theorem 4.4 of [START_REF] Benzoni-Gavage | Multidimensional hyperbolic partial differential equations[END_REF]). At last the "weak=strong" lemma does not require a lot of comments because in the strip geometry there is only one normal variable, namely x d . So that the classical regularization by a mollifier in the tangential variables applies exactly as in the half-space geometry (see [START_REF] Benzoni-Gavage | Multidimensional hyperbolic partial differential equations[END_REF] Section 4.5).

This ends the proof of Theorem 2.1.

4 Saturation of (5), proof of Theorem 2.2

In this section we give a proof of Theorem 2.2. The sketch of proof that we follow here relies on the construction of a geometric optics expansion for (1) and consequently it is probably not the straighest possible proof7 .

However we believe that this sketch of proof is interesting for several reasons. The first one is that giving an approximate solution to (1) is an interesting result in its own. The second is that the construction in itself explains why no new loss of derivative occurs.

And finally, because the construction of the geometric optics expansions for an elliptic frequency (on which the uniform Kreiss-Lopatinskii degenerates) is simpler than the one for a hyperbolic frequency (see [Benc]).

Construction of a geometric optics expansion

Let ζ := (iτ , η) ∈ Υ in all the following we consider the phase function ψ := ψ(t, x ) = τ t + η • x and we consider the following highly oscillating boundary value problem: for 0 < ε 1

           L(∂)u ε = f ε , in Ω, B 0 u ε |x d =0 = g ε 0 for (t, x ) ∈ ω, B 1 u ε |x d =1 = g ε 1 for (t, x ) ∈ ω, u ε |t≤0 = 0 on Γ, (13)
where the source terms in (13) are given by:

f ε (t, x) := e i ε ψ f ev t, x; x d ε + f ex t, x; x d -1 ε , g ε 0 (t, x ) := εe i ε ψ g 0 (t, x ) and g ε 1 (t, x ) := e i ε ψ g 1 (t, x ).
the boundary amplitudes g 0 , g 1 ∈ H ∞ (ω) vanish for negative times while the interior amplitudes f ev ∈ P ev (Ω), f ex ∈ P ex (Ω); the profile spaces defined as follows:

Definition 4.1 (Profile space) We denote by P ev (Ω) (resp. P ex (Ω)) the set of evanescent (resp. explosive) profiles. More precisely, P ev (Ω) (resp.

P ex (Ω)) is the set of functions U = U (t, x; X) ∈ H ∞ (Ω × R + ) (resp. H ∞ (Ω × R -))
vanishing on negative times and for which there exists δ > 0 such that e δX U (t, x; X) ∈

H ∞ (Ω × R + ) (resp. H ∞ (Ω × R -)).
For later use we also introduce the following vectors: We postulate for ansatz:

u ε (t, x) ∼ n≥0 ε n e i ε ψ U ev,n t, x; x d ε + U ex,n t, x; x d -1 ε , (14)
where for all n ≥ 0, U n,ev ∈ P ev and U n,ex ∈ P ex .

Injecting the ansatz (14) in the evolution equation of ( 13) gives (because of the definition of the profile spaces) the cascade of equations:

         L(∂ X )U 0,ev = 0 on Ω × R + , L(∂ X )U 0,ex = 0 on Ω × R -, L(∂ X )U n+1,ev + L(∂)U n,ev = δ n,0 f ev on Ω × R + and ∀n ≥ 0, L(∂ X )U n+1,ex + L(∂)U n,ex = δ n,0 f ex on Ω × R -and ∀n ≥ 0, (15)
where δ •,• stands for the Kronecker symbol and where the operator of derivation in the fast variable L(∂ X ) is given by

L(∂ X ) := A d ∂ X -A (ζ) ,
in which we recall that A (ζ) stands for the resolvent matrix defined in (3).

In order to determine 0 we consider the amplitude of order one. Because of Lemma 4.1 applied to the second equation of (16) (for n = 1) we have (I -P ev )U 1,ev = Q ev (f ev -L(∂)U ev,0 ). Injecting this expression in the boundary condition of (16) (for n = 1) then gives b

• B 0 (P ev U 1,ev ) |x d =X=0 =0 -b • B 0 ∞ 0 e -sA (ζ) Π u (ζ)A -1 d (f ev -L(∂)U 0,ev ) (t, x , 0; s) ds = b • g 0 .
Expanding the left hand side of the previous equation shows that 0 satisfies the transport equation:

b • B 0 I 0 ∂ t + d-1 j=1 b • B 0 I j ∂ j 0 = g 0 in ω, 0 |t≤0 = 0, (18)
where

g 0 := b • g 0 + B 0 ∞ 0 e -sA (ζ) Π u (ζ)A -1
d f ev ds and for all j ∈ 0, d -1 , ζ) e ds (with the convention A 0 = I) .

I j := ∞ 0 e -sA (ζ) Π u (ζ)A -1 d A j e sA (
Reiterating the same proof as the one of [START_REF] Marcou | Rigorous weakly nonlinear geometric optics for surface waves[END_REF], one can show that for all j ∈ 0, d -1 ,

I j = Π u (ζ)∂ ηj Π s (ζ)e = ∂ ηj e(ζ)
(with the convention η 0 = τ ). A direct consequence of Assumption 2.3 is that the I j can be expressed in terms of the (real valued) function θ. More precisely, we obtain that (18) is equivalent to the following transport equation (with real coefficients) A simple integration along the characteristics then determines the value of 0 and the value of the "double" trace U 0,ev (t, x , 0; 0). To complete the determination of U 0,ev it is then sufficient to extend this "double" trace in a single one by setting, for example, U 0,ev (t, x; X) := χ(x d ) 0 (t, x )e XA (ζ) e,

b • β(ζ) ∂ τ θ(ζ)∂ t + ∇ η θ(ζ) • ∇ x 0 = -i g 0 in ω, 0 |t≤0 = 0, (19)
where χ ∈ D(]-1/2, 1/2[ ; R) satisfies χ(0) = 1.

Reiterating essentially the same kind of computations, we can easily construct the higher order amplitudes, namely the (U n,ev ) n≥1 and the (U n,ex ) n≥1 and show that the truncated ansatz ( 14) is a good approximation of the exact solution u ε to (13). More precisely we have the following result: Theorem 4.1 Under Assumptions 2.1-2.2 and 2.3, there exists (U n,ev ) n≥0 ∈ P N ev (resp. (U n,ex ) n≥0 ∈ P N ex ) satisfying the cascade of equations (16) (resp. (17)). Moreover let u ε be the exact solution to (16) and for N 0 ∈ N let

u ε N0 := N0 n=0
ε n e i ε ψ U n,ev t, x;

x d ε + U n,ex t, x; x d -1 ε , (21)
then for all T > 0, α ∈ N d+1 and ε ∈ ]0, 1] we have the estimate

ε |α| u ε -u ε N0 L 2 (Ω T ) α ε N0+3/2 .

Proof of Theorem 2.2

With Theorem 4.1 in hand, the proof of Theorem 2.2 is quite immediate. By contradiction we assume that s = 1 2 -δ for some δ > 0. In (13) let the source terms be given by f ε ≡ 0, g ε 1 ≡ 0 and g ε 0 (t, x ) := εe i ε ψ(t,x ) κ(t, x )b, where κ ∈ D(ω; R). In the following, u ε stands for the exact solution to (13).

By standard interpolation it turns out g ε 0 is O(ε 1/2+δ ) in H s (ω T ). So the energy estimate (7) implies that u ε is O(ε 1/2+δ ) in L 2 (Ω T ). Theorem 4.1 applied to N 0 = 1 combined with the triangle inequality show that u ε 0 (see ( 21)) is also O(ε 1/2+δ ).

Using a technical lemma from [START_REF] Benoit | Geometric optics expansions for linear hyperbolic boundary value problems and optimality of energy estimates for surface waves[END_REF] (more precisely Lemma 4.1) it implies that U 0,ev is zero in L 2 (Ω T ). This is a contradiction with the explicit formula (20) in which 0 can not be zero because it is the solution to the transport equation ( 19) with interior source term g 0 = κ|b| 2 ≡ 0.

  The boundary conditions B 0 and B 1 satisfy the following 1. The boundary condition on {x d = 1} satisfies the uniform Kreiss-Lopatinskii condition that is to say that ∀ζ ∈ Ξ, kerB 1 ∩ E u (ζ) = {0} .2. The boundary condition on {x d = 0} violates the uniform Kreiss-Lopatinskii condition at order one in the elliptic area 4 . More precisely if we denote by Υ := {ζ ∈ Ξ \ ker B 0 ∩ E s (ζ) = {0}} then (a) the weak Kreiss-Lopatinskii condition is satisfied that is ∀ζ ∈ Ξ \ Ξ 0 , ker B 0 ∩ E s (ζ) = {0} , and consequently Υ ⊂ Ξ 0 . (b) We have Υ = ∅ and Υ ⊂ E . (c) Let ζ ∈ Υ then there exist a neighbourhood V of ζ ∈ Ξ, a basis (E 1 , ..., E p )(ζ) of E s (ζ) regular on V , a matrix P = P (ζ) ∈ GL p (C) regular on V and a real valued function θ regular on V such that ∀ζ ∈ V , B 0 [E 1 , ..., E p ] (ζ) = P (ζ)diag(γ + iθ(ζ), 1, ..., 1).

Definition 4. 2

 2 Under Assumptions 2.1, 2.2 and 2.3, let ζ ∈ Υ we can define: • a vector e := e(ζ) ∈ C N \ {0} such that ker B 0 ∩ E s (ζ) = vect {e}. • A vector b := b(ζ) ∈ C p \ {0} such that for all v ∈ E s (ζ) we have b • B 0 v = 0. • Finally a supplementary space Ěs := Ěs (ζ) to vect{e} in E s (ζ) such that B is an isomorphism from Ěs to b ⊥ .

  whereβ(ζ) ∈ C p the first column of P (ζ) (see Assumption 2.3) satisfies b • β(ζ) = 0.

All the results of the paper hold mutatis mutandis if each boundary condition is responsible of a degeneracy of the uniform Kreiss-Lopatinskii condition on an elliptic frequency.

It is in fact zero when the boundary conditions responsible of the failure of the uniform Kreiss-Lopatinskii condition are homogeneous.

Note that the existence of such a configuration implies that N = 2p.

The definition given here differs a little in its formulation (but is equivalent) from the one given in [[ST88] condition (LC)]. The advantage of this formulation will be made precise in Section 4.

Theorem 2.1 and its proof can, under the suitable assumptions, probably be extended to time/space depending coefficients matrices A j , B 0 and B 1 .

We do not expose here the precise functional or the spaces used. However they are exactly the same as the ones used in [[ST88]-Section 1.3] up to the fact that the extra boundary condition on {x d = 1} has to be included. However this boundary condition is harmless because it satisfies the uniform Kreiss-Lopatinskii condition so the standard duality method applies directly.

Indeed the fact that the estimate of[START_REF] Sablé-Tougeron | Existence pour un problème de l'élastodynamique Neumann non linéaire en dimension 2[END_REF] is sharp in the half-space geometry combined with finite speed of propagation and localization arguments make the sharpness of (5) almost readable.

Moreover this construction is more or less contained (in a simpler setting) in the construction of the (Uev,n) n∈N bellow.

Then plugging the ansatz (14) in the boundary conditions of (13) leads to: ∀n ≥ 0 B 0 U ev,n (t, x , 0; 0) + B 0 U ex,n t, x , 0; -1 ε = δ n,1 g 0 on ω, B 1 U ev,n t, x , 1; 1 ε + B 1 U ex,n (t, x , 1; 0) = δ n,0 g 1 on ω.

However by definition of the profile spaces the terms U ex,n t, x , 0; -1 ε and U ev,n t, x , 1; 1 ε are O(ε ∞ ) so that they can be neglected in the previous boundary conditions.

Consequently in the particular setting where the frequency is elliptic we obtain a total decoupling of the cascade of equations. It turns out that the cascades of equations to solve are:

and

on ω and ∀n ≥ 0, U ex,n |t≤0 = 0 on Γ and ∀n ≥ 0.

(17)

In the following the construction of a sequence (U ex,n ) n∈N satisfying (17) will not be described because it is a slight adaptation of the works of [Mar10]- [START_REF] Lescarret | Wave transmission in dispersive media[END_REF] (see also [START_REF] Williams | Nonlinear geometric optics for hyperbolic boundary problems[END_REF]- [START_REF] Williams | Boundary layers and glancing blow-up in nonlinear geometric optics[END_REF] for a construction in which the elliptic modes are treated mode by mode, and not by a one block argument) 8 .

The construction of a sequence (U ev,n ) n∈N satisfying ( 16) is based on the work of [START_REF] Marcou | Rigorous weakly nonlinear geometric optics for surface waves[END_REF] (up to the little modification induced by a non-zero source term on the boundary). However, because this mode carries the loss of regularity of the boundary value problem (1) we sketch in the following the main steps in the construction of U 0,ev .

The key ingredient in the construction of U 0,ev is the following lemma (a consequence of Duhamel formula)):

Then for all F ∈ P ev the equation L(∂ X )U = F on Ω × R + admits a solution U ∈ P ev written under the form U = P ev U + Q ev F .

From Lemma 4.1 we obtain the "polarization" relation on the leading order amplitude that is U 0,ev = P ev U 0,ev . Consequently from the definition of P ev we deduce that to determine U 0,ev it is sufficient to determine U 0,ev (t, x; 0) and that U 0,ev ∈ E s (ζ).

To determine U 0,ev (t, x; 0), we first determine the "double" trace U 0,ev (t, x , 0; 0) and then to extend this trace to x d > 0 we consider it as a boundary layer.

The boundary condition of (17) written for n = 0 implies that U 0,ev (t, x , 0; 0) ∈ ker B 0 so that there exists a scalar function 0 defined on ω such that U 0,ev (t, x , 0; 0) = 0 (t, x )e and the construction of U 0,ev only requires the determination of 0 .


Sneha Singh

Yann Capdeville

Heiner Igel

Correcting wavefield gradients for the effects of local small-scale heterogeneities

Keywords: Site effects, Theoretical Seismology, Rotational Seismology

Local small-scale heterogeneities, like geological inhomogeneities, surface topographies and cavities are known to affect seismic wavefield gradients. This effect is in fact measurable with current instruments. For example, the agreement between data and synthetics computed in a tomographic model is often not as good for rotation as it is for displacement. Here, we use the theory of homogenization to explain why small-scale heterogeneities strongly affect wavefield gradients, but not the wavefield itself. We show that at any receiver measuring wavefield gradient, small-scale heterogeneities cause the wavefield gradient to couple with strain through a coupling tensor J. Furthermore, we show that this J is (1) independent of source, (2) independent of time, but (3) only dependent on the receiver location. Consequently, we can invert for J based on an effective model for which synthetics fit displacement data reasonably well. Once inverted, J can be used to correct rotations at that receiver for other sources. Results of the correction are shown for synthetic rotational receivers and for the ring laser located in Wettzell, Germany. We find that compared to synthetics, the synthetics corrected through J are a better fit to the rotation data. Although results here are derived for rotations, they can be extended to receivers measuring any wavefield gradient.

INTRODUCTION

Seismic sources produce not only translations but also wavefield gradients, i.e. strains and rotations. Therefore, to assess any ground motion completely, the six components of strain and the three components of rotation are as important as the three components of translation [START_REF] Aki | Quantitative Seismology[END_REF]. Collocated measurements of motions can be of great benefit for seismology because they allow us to retrieve instantaneous slowness without using travel times [START_REF] Edme | Local dispersion curve estimation from seismic ambient noise using spatial gradients[END_REF][START_REF] Wassermann | Toward a single-station approach for microzonation: Using vertical rotation rate to estimate Love-wave dispersion curves and direction finding[END_REF][START_REF] Bernauer | Measurements of translation, rotation and strain: new approaches to seismic processing and inversion[END_REF][START_REF] Fichtner | Sensitivity densities for rotational ground-motion measurements[END_REF][START_REF] Igel | Broad-band observations of earthquake-induced rotational ground motions[END_REF][START_REF] Igel | Rotational motions induced by the M8.1 Tokachi-oki earthquake[END_REF][START_REF] Mikumo | Determination of local phase velocity by intercomparison of seismograms from strain and pendulum instruments[END_REF], to constrain seismic source properties [START_REF] Donner | Retrieval of the seismic moment tensor from joint measurements of translational and rotational ground motions: Sparse networks and single stations[END_REF][START_REF] Donner | Inversion for seismic moment tensors combining translational and rotational ground motions[END_REF][START_REF] Reinwald | Improved finite-source inversion through joint measurements of rotational and translational ground motions: a numerical study[END_REF][START_REF] Bernauer | Reducing nonuniqueness in finite source inversion using rotational ground motions[END_REF][START_REF] Takeo | What can be learned from rotational motions excited by earthquakes?[END_REF] and to estimate parameters like propagation direction [START_REF] Igel | Broad-band observations of earthquake-induced rotational ground motions[END_REF][START_REF] Cochard | Rotational motions in seismology: theory, observation, simulation[END_REF]) and rupture history [START_REF] Takeo | What can be learned from rotational motions excited by earthquakes?[END_REF]. Concurrent measurements of all three kinds of motions also help us remove the coupling that one motion introduces on the other, for example, horizontal rotations on translation (Aki & Richards 2002, p. 604) and vertical rotation on translation [START_REF] Graizer | The response to complex ground motions of seismometers with Galperin sensor configuration[END_REF][START_REF] Trifunac | A note on the useable dynamic range of accelerographs recording translation[END_REF]. Further applications of collocated measurements have been discussed for example, in [START_REF] Schmelzbach | Advances in 6C seismology: Applications of combined translational and rotational motion measurements in global and exploration seismology[END_REF], [START_REF] Li | Tutorial on rotational seismology and its applications in exploration geophysics[END_REF] and [START_REF] Igel | Rotational seismology: Theory, instrumentation, observations, applications[END_REF].

Wavefield gradients themselves are important to deduce the stresses caused by seismic waves [START_REF] Spudich | Transient stresses at Parkfield, California, produced by the M7.4 Landers earthquake of June 28, 1992: Observations from the UPSAR dense seismograph array[END_REF], to enable S wave interpretation [START_REF] Sollberger | Singlecomponent elastic wavefield separation at the free surface using source-and receiver-side gradients[END_REF] and in earthquake engineering to assess seismic hazards [START_REF] Basu | Characterizing rotational components of earthquake ground motion using a surface distribution method and response of sample structures[END_REF][START_REF] Castellani | On the rotational components of seismic motion[END_REF][START_REF] Abdel-Ghaffar | Torsional earthquake response of suspension bridges[END_REF][START_REF] Trifunac | A note on rotational components of earthquake motions on ground surface for incident body waves[END_REF]. Despite all these applications, measurements of wavefield gradients are relatively few. They were primarily estimated using dense arrays [START_REF] Huang | Ground rotational motions of the 1999 Chi-Chi, Taiwan earthquake as inferred from dense array observations[END_REF][START_REF] Bodin | Dynamic deformations of shallow sediments in the valley of Mexico, part I: Three-dimensional strains and rotations recorded on a seismic array[END_REF][START_REF] Spudich | Transient stresses at Parkfield, California, produced by the M7.4 Landers earthquake of June 28, 1992: Observations from the UPSAR dense seismograph array[END_REF] for a very long time because seismic instruments to measure them were either not reliable enough or sensitive enough in the required frequency bands, thus leading (Aki & Richards 1980, p. 489) to note that seismology still awaits a suitable instrument for such measurements. Nowadays, rotations can be measured using liquid-based motion sensors [START_REF] Egorov | Self-noise of the MET angular motion seismic sensors[END_REF][START_REF] Huang | Molecular electric transducers as motion sensors: A review[END_REF], magneto-hydrodynamic sensors [START_REF] Pierson | Advances in rotational seismic measurements[END_REF], ring laser gyroscopes [START_REF] Schreiber | Ring laser measurements of ground rotations for seismology[END_REF] Correcting wavefield gradients for the effects of local small-scale heterogeneities 3 [START_REF] Pancha | Ring laser detection of rotations from teleseismic waves[END_REF][START_REF] Mcleod | Comparison of standard and ring laser rotational seismograms[END_REF], mechanical sensors [START_REF] Brokešová | Rotaphone, a self-calibrated six-degree-of-freedom seismic sensor and its strong-motion records[END_REF]) and more recently, using fiber-optic gyroscopes [START_REF] Bernauer | BlueSeis3A: Full characterization of a 3C broadband rotational seismometer[END_REF][START_REF] Lindner | Seafloor ground rotation observations: Potential for improving signal-to-noise ratio on horizontal OBS components[END_REF][START_REF] Kurzych | Fibre optic system for monitoring rotational seismic phenomena[END_REF]. Similarly, strains can be measured using Distributed Acoustic Sensing (DAS) [START_REF] Willis | Quantitative quality of distributed acoustic sensing vertical seismic profile data[END_REF] and Distributed Vibration Sensing (DVS) [START_REF] Dean | The effect of gauge length on axially incident P-waves measured using fibre optic distributed vibration sensing[END_REF].

Although considerable effort has been put into devising more robust and portable instruments, little attention has been paid to yet another issue concerning wavefield gradients. It has been shown that local small-scale heterogeneities, such as cavities and discontinuities, topographies as well as geological heterogeneities strongly affect strain [START_REF] Gomberg | The accuracy of seismic estimates of dynamic strains: An evaluation using strainmeter and seismometer data from Piñon Flat Observatory, California[END_REF][START_REF] Berger | An analysis of tidal strain observations from the United States of America II. The inhomogeneous tide[END_REF]) and rotation [START_REF] Van Driel | Strain rotation coupling and its implications on the measurement of rotational ground motions[END_REF][START_REF] Kohl | Measurement and interpretation of tidal tilts in a small array[END_REF][START_REF] Harrison | Cavity and topographic effects in tilt and strain measurement[END_REF][START_REF] King | Tidal tilt measurement in europe[END_REF]. This effect of heterogeneities on wavefield gradient is measurable with our current sensors. This leads to a number of consequences, such as 1) disagreement between wavefield gradient data and synthetics computed in a tomographic model, 2) dissimilarity between wavefield gradients measured at a point and wavefield gradients estimated at the point using an array, and 3) discrepancy between data observed among closely-spaced instruments. Such consequences may lead to an erroneous interpretation of the data.

The purpose of this paper is twofold. First, we study why wavefield gradients are affected by local small-scale heterogeneities. This will help us understand how a wavefield and its observables really interact with heterogeneities on all scales, which will eventually help us in better interpreting seismic data. Second, we investigate if it is possible to correct wavefield gradients for this effect. This will ensure a better agreement between wavefield gradient data and the synthetics that have been corrected. This will also help align the wavefield gradients measured at a point and wavefield gradients estimated at the point using an array.

The paper is organised as follows: we first describe an example in section 2 that will enable us to generalize the effect of small-scale heterogeneities on wavefield gradients. Using rotations and strains as our measurements, we show that indeed, small-scale heterogeneities strongly affect wavefield gradients. Then we approach the physics of this problem following the theory of homogenization. We demonstrate that there is a possibility of correcting for this effect. Finally, in section 3, we show the results of this correction at a few synthetic rotational receivers before moving on to the G Ring laser in section 4. Although examples shown in 3 and 4 concern rotation, the results can be extended to any other wavefield gradient and their time derivatives.

EFFECTS OF SMALL-SCALE HETEROGENEITIES ON WAVEFIELD GRADIENT

MEASUREMENTS

In this section, we begin with a simple numerical example where we show that the effect of smallscale heterogeneities is of a different nature on wavefield gradients, such as rotation and strain compared to the wavefield, such as displacement, velocity, and acceleration. This will be followed by theory in section 2.2 where we make use of homogenization to explain the observed effect and by a discussion in section 2.3 where we demonstrate the correction of this effect.

Numerical Example

We consider a 2D elastic case where the computational domains are two 20×20 km 2 squares as shown in Figure 1. These squares, except for spatially varying small-scale heterogeneity at the centre of the second (Fig. 1; model on the top-right panel), are otherwise homogeneous. The P and S velocities in the homogeneous area of the squares are 5 kms -1 and 3.2 kms -1 respectively.

The spatial extension of the heterogeneity at the center of the second square is 200×200 m 2 . Its geometry is kept simple in order to make both meshing and waveform modeling easier. The velocities in each element of the heterogeneity differ by ±50% from the surroundings. This velocity contrast is strong but realistic for such a scale range. For example, the velocity contrast between granite and unconsolidated sediments is usually even larger. Finally, four receivers are placed in both squares such that for the heterogeneous square, two receivers lie just inside the heterogeneity and the other two just outside (Fig. 1; bottom panel). The distance between the receivers is around 15 m. All boundaries are absorbing.

In both cases, we compute the wavefield induced by a horizontal double-couple source. The Bottom: V p , V s and ρ at the heterogeneity. The length of the minimum wavelength for the simulation is shown in the leftmost panel.

than λ m . Using a 2D Spectral Element [START_REF] Komatitsch | The spectral-element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF] program, displacements, velocities, accelerations, shear strain and rotation measured by the four receivers in the two models are computed. These are shown in Figure 2. Because the mesh used for the computation honors the interfaces of the heterogeneity, the solutions thus obtained are considered to be accurate.

It follows from Figure 2 that while all four receivers measure almost the same displacement, velocity and acceleration in both models, this is not true for rotation and strain. The four receivers measure a very similar rotation and strain in the homogeneous case, but in the presence of a heterogeneity, though much smaller than the minimum wavelength, the measurements are substantially different. The difference can be observed in both amplitudes and phases, and it depends on the rel- 5.0 x 10 -3 9.0 x 10 -3 3.0 x 10 -3 11.0 x 10 -3 Correcting wavefield gradients for the effects of local small-scale heterogeneities 7 ative location of the receiver to the heterogeneity. Because our current sensors are sensitive enough to measure this difference, interpretations derived from the four receivers would be incoherent. For example, if these rotations were to be used to determine instantaneous slowness, one would arrive at four different values of slowness. Note that this result is not specific to just the quantities shown.

The conclusion would be the same for any measurements of wavefield (and its time derivatives, such as displacement, velocity and acceleration) and for any measurements of wavefield gradient (and its time derivatives, such as strain, strain rate, rotation and rotation rate). It is also worth noting that acceleration is insensitive to the heterogeneity while rotation rate is. Because rotation rates are related to local phase velocities through transverse accelerations [START_REF] Cochard | Rotational motions in seismology: theory, observation, simulation[END_REF][START_REF] Igel | Rotational motions induced by the M8.1 Tokachi-oki earthquake[END_REF][START_REF] Pancha | Ring laser detection of rotations from teleseismic waves[END_REF], this might seem inexplicable. However, the relation is only valid in the absence of heterogeneities, and as Figure 2 shows, it no longer holds when the receiver is placed at a heterogeneity. Similarly, Brokešová & Málek (2015a) and Brokešová & Málek (2015b) have shown that under certain conditions, velocities might be as important as accelerations for relating rotations to translations. However, as Figure 2 shows, the effect on rotations does not originate from velocities either. From this simple numerical example, it can be concluded that small-scale heterogeneities have almost no effect on wavefield measurements but have a much stronger effect on wavefield gradient measurements. This observation is an issue for any wavefield gradient based instrument as the measurements can be tainted by nearby heterogeneities. A more intensive set of experiments would have shown that the amplitude of this effect depends mainly on the velocity contrast and not the heterogeneity's spatial extent with respect to the wavelength. In the following section, we examine this effect theoretically.

Theory

Homogenization deals with obtaining macroscopic, or homogenized or effective equations for systems with microscopic structures. In seismology, a microscopic medium is characterized by variations in elastic properties at scales much smaller than the minimum wavelength in the medium. By using homogenization, all these microscopic structures are smoothened out, i.e. the medium with rapidly-varying structures can be substituted with another where the structures are only slowlyvarying. To this end, homogenization has been beneficial in forward problem to study wave propagation in a fine-layered media [START_REF] Backus | Long-wave elastic anisotropy produced by horizontal layering[END_REF], to allow sparser meshing to gain computational speed [START_REF] Cupillard | Non-periodic homogenization of 3-D elastic media for the seismic wave equation[END_REF]Capdeville et al. 2010a) as well as in the inverse problem to help interpret the results of full waveform inversion [START_REF] Capdeville | Elastic full waveform inversion based on the homogenization method: Theoretical framework and 2-D numerical illustrations[END_REF]. Because homogenization allows us to understand how different scales interact with the wavefield, it identifies itself as a natural choice to study the effects of small-scale structures on the gradients of displacement.

In the following, we give an intuitive construction of the homogenization problem. For simplicity, we consider a 2D periodic media where the properties repeat after certain intervals. For more extensive description on homogenization, Capdeville et al. (2010b,a); [START_REF] Cupillard | Non-periodic homogenization of 3-D elastic media for the seismic wave equation[END_REF] are recommended for 1D, 2D and 3D respectively.

We consider a wavefield u(x, t) which is the solution of elastic wave equation for an external source f (x, t) with a maximum frequency f max . Because most geophysical media have no specific scale separation, we manually assign a separation scale λ. All spatial variations of the wavefield at scales below λ are considered small-scales whereas all those above λ are considered large-scales.

Accordingly, we introduce the parameter ε,

ε = λ λ m , (1)
where λ m denotes the minimum spatial wavelength of the wavefield u.

To separate large scales from small scales explicitly, we introduce a new space variable y,

y = x ε . ( 2
)
y represents all the small-scales and x represents all the large-scales. It is then assumed that x and y can be treated as independent variables. This two-scale assumption can be mathematically justified by demonstrating its convergence towards the solution of a one-scale variable (for example, Dal Maso 1993; Allaire 1992). Consequently, all gradients with respect to x are redefined

Correcting wavefield gradients for the effects of local small-scale heterogeneities 9 as

∇ x → ∇ x + 1 ε ∇ y , (3)
where

∇ x = t (∂ x1 , ∂ x2 ) and ∇ y = t (∂ y1 , ∂ y2 ).
The aim of homogenization is to find solutions to the wave equations as an asymptotic ex-

pansion in ε [see Appendix A].
As such, the equations for true displacement u ε and true gradient (∇u) ε take the form

u ε (x, y, t) = u 0 (x, y, t) + εu 1 (x, y, t) + ε 2 u 2 (x, y, t) + ... (∇u) ε (x, y, t) = (∇u) 0 (x, y, t) + ε(∇u) 1 (x, y, t) + ε 2 (∇u) 2 (x, y, t) + ... (4)
In ( 4), the superscripts are powers in ε but not in the u and (∇u).

Solving the series of homogenization equations [see Appendix A], it can be shown that,

u 0 (x, y, t) = u 0 (x, t), (5)
u 1 (x, y, t) = χ(y) : (u 0 (x, t)). (6)
In ( 6), denotes the strain ([

(u 0 )] ij = 1 2 (∂ i u 0 j + ∂ j u 0 i ))
and is not to be confused with the homogenization parameter ε from (1). χ, the first-order corrector, is a third-order tensor. It is a solution of the cell equations [see Appendix A]. Note that χ does not depend on time or sources.

Therefore, to leading order, the true displacement u ε can be rewritten as

u ε (x, y, t) = u 0 (x, t). (7)
The true gradient (∇u) ε to the leading order is then

(∇u) ε (x, y, t) = (∇ x u 0 )(x, t) + ∇ y χ(y) : (u 0 (x, t)). (8)
This implies that for rotation in 2D, which is a scalar,

ω ε = θ ij (∂ i u j ) 0 = ω 0 + J ij : ij (u 0 ), (9)
where J ij = θ mn ∂ ym χ ijn and θ represents the Levi-Civita symbol. Similarly, for rotation in 3D,

ω ε i = θ ijk (∂ j u k ) 0 = ω 0 i + J ijk : jk (u 0 ), (10)
where J ijk = θ imn ∂ ym χ njk and θ represents the Levi-Civita symbol. And for a rotational receiver on the surface of a 3D domain under free surface conditions,

ω ε i = ω 0 i + Jijk : σ jk , (11)
where Jijk = J imn s mnjk , σ is the stress tensor and s is the inverse of the elastic tensor c.

In (4), u ε is the true displacement and u 0 is the zeroth order effective displacement. Similarly, (∇u) ε represents the true gradient of displacement and (∇u) 0 represents the zeroth order effective gradient of displacement. True quantities correspond to real data that can be collected in the field. In theory, it is possible to model them using wave equations but because they depend on all possible scales, they require an infinite scale knowledge of the elastic properties of the Earth. This is not feasible in practice. On the other hand, for a given frequency band, effective quantities can be modeled using the same equations but using homogenized elastic properties which do not contain any small scales. Thus, effective quantities can be described with a finite set of parameters.

Retrieving these parameters is what is actually done in full waveform inversion and tomography imaging [START_REF] Capdeville | Elastic full waveform inversion based on the homogenization method: Theoretical framework and 2-D numerical illustrations[END_REF]. Most of the tomographic models that exist are therefore actually an approximation to the homogenized Earth.

Comparing the terms in (7), we note that the terms on the right depends only on x, i.e. largescales. This means that effective displacement can be used as good approximation for true displacement without worrying about small-scales. This is a well known result from Chapman (2004, p. 277) and [START_REF] Backus | Long-wave elastic anisotropy produced by horizontal layering[END_REF]. One consequence of this result is that the agreement between data and synthetics is reasonably good for displacement, velocity and acceleration. This is also consistent with the results obtained from section 2.1 and with many examples concerning real data. For instance, very long period surface waves can be modelled using PREM with good accuracy even though the crust is highly heterogeneous. The large amplitude and small-scale heterogeneity of

Correcting wavefield gradients for the effects of local small-scale heterogeneities 11 the crust is inconspicuous to the displacement field, and one does not need complete information on the crust for its modeling.

On the contrary, this is not true for any gradient of displacement. The terms in (8) depend on both x and y, i.e. both large-scales and small-scales. Due to small-scales, the true gradient of displacement now couples with strain of effective displacement. The coupling is due to the firstorder corrector χ, therefore, we refer to this coupling as 'correction'. Similarly, we refer to J as the coupling tensor.

The appearance of correction in (8) implies that setting up an inversion problem with wavefield gradient is impossible because observations rely on all scales. Every detail in all scales would be necessary to set up the inversion problem. For the case of displacement, velocity or acceleration where effective quantities are a good approximation to the true quantities, this is not the case and setting up an inversion problem is not an issue. For the wavefield gradients, however, correction needs to be accounted for and this is the purpose of the next section.

Inversion of the corrector

We know from ( 9), ( 10) and ( 11) that the true wavefield gradient is coupled with strain of the effective displacement through the tensor J. Moreover, J depends neither on time nor on source. It depends only on the location of the receiver. If the elastic properties at the receiver are known well enough, J can be computed by solving the cell problem (Capdeville et al. 2010a,b). This is not the case in practice. Nevertheless, we can invert for J provided that we have an effective model that fits the displacement (or velocity or acceleration) data reasonably well. The obvious choice for this effective model is the homogenized model because this is what the waves can see [START_REF] Capdeville | Elastic full waveform inversion based on the homogenization method: Theoretical framework and 2-D numerical illustrations[END_REF]. This effective model is usually not known, but any tomographic model can be a good approximation as we shall see from the synthetic examples in section 3. For real data, any

Earth model which fits the displacement (or velocity or acceleration) data reasonably well can be used. This is similar to an earthquake moment tensor inversion where an Earth model is necessary to compute Green's functions.

Because all examples in this paper are carried out for rotation, we start from the equation

ω ε (x r , t; x s ) = ω 0 (x r , t; x s ) + J(y r ) : 0 (x r , t; x s ), (12)
where ω ε is the rotation data, ω 0 and 0 are the rotation and the strain computed in the effective model and J is the unknown. x s denotes the location of the source and y r = xr ε denotes the location of the receiver. We invert for J by solving the over-determined linear problem in the least-squares sense [START_REF] Tarantola | Generalized nonlinear inverse problems solved using the least squares criterion[END_REF]. We seek to minimise the L 2 misfit between the synthetics and the data through

J ∼ ( t G G) -1 G δd, (13)
where δd = ω ε -ω 0 and G is the matrix where all strains of the effective displacement are stacked.

The inversion process is straight-forward: we need only one simulation per source to assemble G. Then, we apply Equation ( 13) with no modification as the inversion is stable and no damping is necessary. We invert for 3 parameters in 2D and for 18 parameters for the 3 components of rotation in 3D. Like the instrument response, the corrector can be used to separate unwanted site effects from data.

Once we invert for J, we need to evaluate the quality of our inversion results. One way to do so would be to compare the J computed through inversion to the one obtained by solving the series of homogenization equations [see the cell problem in Appendix A]. However, this solution is only available when all fine-scales are perfectly known, i.e only for synthetic tests. Therefore, for all our examples, we start by dividing the data into two pools: the inversion data set and the validation data set. We use the inversion data set to perform the inversion and the validation data set to assess the accuracy of J. In all examples that follow, we use the J to correct the synthetic rotations.

Therefore, a better fit between observed data and corrected synthetics, especially for events not used in the inversion, implies that J is a good measure of the coupling due to small-scales at the receiver.

Correcting wavefield gradients for the effects of local small-scale heterogeneities 13 V p , V s and ρ in the domain. The P and S velocities in the outer domain are 5 kms -1 and 3.2 kms -1 , and velocities for each element inside the heterogeneity are ± 50% of the surrounding domain.

SYNTHETIC TESTS

To test inversion of the coupling tensor J, we first conduct synthetic tests. In section 3.1, we invert for J for a rotational receiver in 2D. This will be followed in section 3.2 where we invert for J for a rotational receiver placed on the free surface of a 3D domain. The set up that has been considered in both cases is quite simple in the sense that it does not take into account factors like attenuation and source uncertainties. All simulations are conducted using the Spectral Element Method (SEM) that solves the elastic wave equation.

In both cases, the sources are either forces or moment tensors but we always choose the Ricker wavelet as the source-time function.

Inversion of J for a receiver in a 2D volume

Our first model is a 2D domain whose P and S velocities are 5 kms -1 and 3.2 kms -1 respectively.

The size of the domain is 20×20 km 2 . A heterogeneity of size 200×200 m 2 is at the center of the domain and the P and S velocities in each element in the heterogeneity are ± 50% of the surrounding domain. The geometric configuration of the model is shown in Figure 3.

Because the rotational receiver has been placed in the heterogeneity, we expect that the heterogeneity affects the rotation measurements. We therefore invert for the coupling tensor J. To this end, we first perform three simulations in this medium with sources at positions shown by the stars a, b and c respectively in Figure 3 to generate the data (Fig. 4; black lines). These rotations make up the term on the left in Equation 12. The central frequency of all three sources is 1.5 Hz (maximum frequency 4.1 Hz). Two of the sources considered are double-couple, while one is a force. Although the heterogeneity is much smaller than λ m , the spectral element mesh for the simulation honors all interfaces of the heterogeneity. The rotations computed are therefore considered accurate.

Next, we need to choose an effective model for the inversion. Certainly, we could use the homogenized version of the original medium but with an underlying assumption that the heterogeneous area has negligible effect on the effective model, we use a homogeneous model instead. This assumption is valid because the displacements in both media are in good agreement [see Appendix C]. Moreover, it is close to realistic situations where tomographic models are far from perfect.

We use the same three sources as before in this homogeneous model and compute the rotations and strains at the receiver. These are the two known quantities on the right in Equation 12. The rotations measured in the homogeneous model are shown in Figure 4 (green lines). Comparison of the rotations measured in the heterogeneous and homogeneous case in the figure reveals a significant difference in all four panels. Clearly, the difference is due to the presence of the heterogeneity and it depends on the path that the waveform takes through the heterogeneity.

We use the data and synthetics thus obtained to invert for J. The fact that we use SEM to both generate data and to invert for J is not an inversion crime here. Indeed, it would be an inversion crime if the data were generated using the theory of homogenization and the inversion performed with the same mesh. Here, we generate the data using a fine mesh with an element size of 100 m, i.e. without the theory of homogenization, but carry out the inversion on a sparse mesh with an element size of 200 m.

Once we obtain J, we use it to correct the synthetics, i.e. we add the last term in Equation 12to the rotations measured in homogeneous medium. These are shown in Figure 4 (red lines). It can be seen that the corrected synthetic rotations are a much better fit to the heterogeneous rotations, but this is unsurprising because these data have been used for the inversion. Therefore, to truly render Heterogeneous and homogeneous rotations are shown in black and green respectively. Rotations that have been corrected using the coupling tensor J are shown in red.

its accuracy, we use this J to a new set of data which has not been used for the inversion (Fig. 4;

Panel d) as validation; the new source positioned at the location shown by the star d in Figure 3.

The new source is a vertical force. As expected, it can be seen that J is able to recover the rotation in the heterogeneous medium.

Inversion of J for a receiver on a 3D surface

Our second model is a 3D domain whose P and S velocities are 5 kms -1 and 3.2 kms -1 respectively. The domain size is 40×40×10 km 3 . A heterogeneity is placed at the center of the surface, with P and S velocities for each element in the heterogeneity chosen randomly and independently but within 25% of the surrounding domain. The spatial extent of the heterogeneity is 1334×1334×222 m 3 . The geometric configuration of the model is shown in Figure 5.

Because the rotational receiver here is surrounded by the heterogeneity, the rotation measurements are most likely affected by the heterogeneity. We therefore seek for the coupling tensor J for Bottom: Zoom in into the heterogeneity from the top showing the V p , V s and ρ in the domain. The P and S velocities in the outer domain are 5 kms -1 and 3.2 kms -1 , and velocities for each element inside the heterogeneity are chosen randomly and independently, but within 25% of the surrounding domain.

this receiver. Here, we only invert for J for one component of the rotation, i.e. rotation around the z-axis. This allows us to invert for just 6 parameters instead of 18. Note that in this case -where the receiver has been placed at the surface -the free surface condition also allows us to further reduce the number of scalars to be inverted for to 3 (Equation 11). To this end, we perform simulations for three different sources at locations shown by the stars a, b and c in Figure 5 in the original medium to generate data (Fig. 6; black lines). These rotations form the term on the left in Equation 12. All sources have a central frequency of 0.5 Hz (maximum frequency 1.4 Hz). Two of the sources are forces, while the third one is a moment tensor. Once again, the mesh honors all boundaries and therefore, the solutions obtained are considered accurate.

Similar to section 3.1, we carry out the inversion using the homogeneous medium as an ap-

Correcting wavefield gradients for the effects of local small-scale heterogeneities 17

proximation for the effective one. We use the same three sources in the homogeneous medium to generate rotations and strains. These are the two known quantities on the right in Equation 12. The rotations in this tomographic medium are shown in Figure 6 (green lines). Comparing the rotations in the heterogeneous and the homogeneous media, we see that the waveform is altered in all four cases owing to the heterogeneity.

We use the data and synthetics thus derived from these simulations to invert for the coupling tensor J which is in turn used to correct the synthetic rotations (Fig. 6; red lines). The corrected synthetic rotations are a much better fit to the heterogeneous rotations, which is expected because these data have been used for the inversion. Finally, we perform one last simulation with the source, which is a moment tensor, at location shown by the star d in Figure 5 to validate the coupling on a new set of data (Fig. 6; Panel d), one that has not been used previously for the inversion. Doing this allows us to estimate the accuracy of J. Unsurprisingly, the recovery of the heterogeneous rotation is achievable through J.

It is important to note that the correction is not always perfect, i. e. the corrected rotations do not always overlap the heterogeneous rotations. This is especially evident in panels a and d in Figure 6. The reasons for this are 1) the homogeneous model used here is not a perfect approximation to the homogenized or effective model and 2) the method, being asymptotic, implies some errors.

REAL DATA

Ever since it's operation in 2001, the G Ring laser in Wettzell, 150 km north-east of Munich (Fig. 8; the seismometer WET and the G Ring laser are collocated), has been successfully recording the vertical component of rotation rates that are generated by seismic sources [START_REF] Suryanto | First comparison of array-derived rotational ground motions with direct Correcting wavefield gradients for the effects of local small-scale heterogeneities 29 ring laser measurements[END_REF][START_REF] Igel | Rotational motions induced by the M8.1 Tokachi-oki earthquake[END_REF][START_REF] Igel | Broad-band observations of earthquake-induced rotational ground motions[END_REF]. Using two counter propagating beams of wavelength 633 nm to generate a beat frequency, it measures rotation rates as small as 10 -10 rads -1 Hz -1/2 [START_REF] Schreiber | The geosensor: A new instrument for seismology[END_REF].

Here, we invert for the coupling tensor J at this ring laser for long period teleseismic events. In section 4.1, we use the normal mode summation program to compute velocities and rotation rates in PREM [START_REF] Dziewoński | Preliminary reference earth model[END_REF]. We show that the agreement between the rotation rates measured at the ring laser and the rotation rates computed in PREM is less compared to the one Heterogeneous and homogeneous rotations are shown in black and green respectively. Rotations that have been corrected using the coupling tensor J are shown in red.

between velocities. We invert for the coupling tensor J and correct the rotation rates computed in PREM to improve the agreement. In section 4.2, we show that due to small-scale heterogeneities at the receiver, observations of rotation rates and array-derived estimations of rotation rates at the ring laser are not the same. We invert for the coupling tensor J and correct the array-derived rotation rates such that they are in better agreement to the observed rotation rates.

We use the definition of L2 norm defined below to quantify the misfit between the observed and synthetic rotation rates.

ξ := tmax 0 ( ω -ωdata ) 2 (t) dt ( ωdata ) 2 (t) dt ; ( 14
)
where ω is the synthetic vertical rotation rate and ωdata is the observed vertical rotation rate. Note that this misfit takes into account all parts of the seismogram.

Correcting wavefield gradients for the effects of local small-scale heterogeneities 19

4.1 Inversion for J at the Ring laser, Wettzell using PREM

We use Obspy [START_REF] Krischer | Obspy : a bridge for seismology into the scientific Python ecosystem[END_REF]) to get rotation rates recorded by the ring laser and velocities recorded by the collocated seismometer WET for all events between January 2007 and June 2018 and within M w 7.0 and 9.0. After download, all data are corrected for instrument response. Next, we compute synthetic velocities and rotation rates by using a normal mode summation program to solve the wave equation in a PREM with ocean and anisotropy. For this, centroid-moment tensor solutions from the GCMT catalog [START_REF] Ekström | The global CMT project 2004-2010: Centroid-moment tensors for 13,017 earthquakes[END_REF]) are used with a step function as the sourcetime function.

We compare the observed and synthetic velocities in order to remove observed data with no visible earthquake signal. Apart from that, no prior selection of events is done. A figure showing the source-receiver distribution in shown in Appendix D. Comparing the fit between observed and synthetic velocities to the fit between observed and synthetic rotation rates for all events within frequencies 0.002-0.009 Hz, we find that the rotation rates are less in agreement compared to the velocities [see Appendix C]. Often, the misfits as defined in ( 14) are about 0.3 between the velocities but are greater than 0.7 between the rotation rates. We suspect this to be due to small-scale heterogeneities at the ring laser and therefore, invert for J using 32 events. The J thus obtained is in turn used to correct the synthetic rotation rates. One of such events is shown in Figure 7 (top panel). A validation data set is then set up to help access the quality of our results. It includes 32 events not used in the inversion. One of the events from this validation data set is shown in Figure 7 (bottom panel). For both panels, we see distinct arrivals of direct Love waves and indirect Love waves. We also see that the corrected rotation rates (red lines) are a better fit to the observed rotation rates (black lines) than the synthetic rotation rates (green lines). Misfits for about 72% of the events in the validation data set improve after the introduction of J, proving that indeed, J is a good indication of the coupling. For example, the misfit for the event shown in the bottom panel of Figure 7 reduced from 0.86 to 0.63. Data, synthetic and corrected rotation rates are given in black, green and red respectively.

4.2 Inversion for J at the Ring laser, Wettzell using array-derived estimates

In the absence of small-scale heterogeneities, wavefield gradients estimated at a point using the finite-difference of array data are the same as wavefield gradients measured at the point. In practice, this is not the case. Wavefield gradients estimated from an array have their coupling tensor although it appears in a slightly different form [see Appendix B] than that of wavefield gradients measured at a point. In the case of rotation, the two coupling tensors are related through the following equations:

ω(x, y, t) = ω 0 (x) + J(y) : (u 0 (x, t)), ω AD (x, y, t) = ω 0 (x) + J AD (y) : (u 0 (x, t)),

where ω 0 is the rotation in an effective model. ω is the rotation measured at a point and J is its coupling tensor. Similarly, ω AD is the rotation estimated at a point using array derived measurements

Correcting wavefield gradients for the effects of local small-scale heterogeneities 21 and J AD is its coupling tensor. These equations lead to ω(x, y, t) = ω AD (x, y, t) + J † (y) : (u 0 (x, t)),

where J † (y) = J(y) -J AD (y).

For our purpose, we use velocities recorded at stations WET, ROTZ and GRC1 (Fig. 8) to estimate array-derived rotation rates at the ring laser. Each station is a broadband STS-2 seismometer that measures velocity with a sampling rate of 20 Hz. The distance between each station is in the order of 100 km. Assuming an average phase velocity of 3 kms -1 , we look at all events between January 2007 and June 2018 for frequencies between 0.0005 and 0.003 Hz following [START_REF] Poppeliers | The effects of measurement uncertainties in seismic-wave gradiometry[END_REF] and [START_REF] Langston | Wave gradiometry in two dimensions[END_REF].

We select only 7 events due to their good signal to noise ratio. Observed, array-derived and corrected rotation rates are given in black, green and red respectively. same procedure as in section 4.1 and use 6 events for the inversion of J † . We then use J † to correct the array-derived rotation rates. One of these events is shown in Figure 9 (top panel). The validation data set in this case consists of 1 event that has not been used for the inversion. This event is shown in Figure 9 (bottom panel). For both events, we see distinct arrivals of direct Love waves and indirect Love waves phases. Furthermore, we see that the corrected rotation rates (red lines) are in better agreement to the observed rotation rates (black lines) than the array-derived rotation rates (green lines). In this case, the misfit for the validation data set was 0.66 originally.

After the introduction of the coupling tensor however, the misfit improved to 0.45.

Note that usually, when estimating array-derived rotations, the effect of small-scale structures is often ignored. In that case, many stations are needed to average out the effect and minimize the coupling tensor J. Here, we explicitly account for the effect and therefore, it is futile to increase the number of stations.

In this study, through a simple example, we showed that local small-scale heterogeneities have little or almost no effect on the wavefield (displacement, velocity and acceleration) but a strong effect on its gradients (strain, strain rate, rotation and rotation rate). This stands in line with previous studies conducted by [START_REF] Gomberg | The accuracy of seismic estimates of dynamic strains: An evaluation using strainmeter and seismometer data from Piñon Flat Observatory, California[END_REF]; [START_REF] Kohl | Measurement and interpretation of tidal tilts in a small array[END_REF]; [START_REF] Lennon | A discussion on the measurement and interpretation of changes of strain in the Earth -Earth tides and their place in geophysics (summary)[END_REF]. Local small-scale heterogeneities can cause both phase and amplitude of the waveform to alter, thus causing a disagreement between the observed data and synthetics. Consequently, these effects are obstacles in direct applications of the measurements, for example, when rotations are used to determine phase velocities [START_REF] Igel | Rotational motions induced by the M8.1 Tokachi-oki earthquake[END_REF] or strains are used to deduce stresses generated by earthquakes [START_REF] Spudich | Transient stresses at Parkfield, California, produced by the M7.4 Landers earthquake of June 28, 1992: Observations from the UPSAR dense seismograph array[END_REF]. Similarly, the free surface condition may not be valid [START_REF] Gomberg | The accuracy of seismic estimates of dynamic strains: An evaluation using strainmeter and seismometer data from Piñon Flat Observatory, California[END_REF].

Through the theory of homogenization, we introduced the coupling tensor J and showed that at any receiver measuring wavefield gradient, small-scale heterogeneities cause the wavefield gradient to couple with the strain of effective displacement through J. Therefore, a knowledge of J is needed to correct for this coupling. One option to determine J would be through direct computation, but this is only feasible if all fine-scales at the receiver are well-known. The second option is to determine J through linear inversion and it is indeed possible to perform this inversion because J does not depend on source or time. It is a characteristic of the receiver location. One necessary condition for this inversion however is that an effective model that fits displacement data well is needed. This is not an issue for most synthetic cases and for low frequency data where PREM can be used as the effective model.

We then inverted for J for rotational receivers in two synthetic cases, one in 2D and one in 3D.

In both these examples, we used the J to correct synthetic rotations and saw that the corrected synthetic rotations were a better match to the synthetic data as compared to synthetic rotations. This was also true for rotations that were not used for inversion. We conclude that although the effect of small-scale heterogeneities can be substantially large, it can be corrected for. The setting for

our synthetic examples was one where the sources are either forces or moment tensors; however, further tests in different settings reveal the same conclusion.

Next, we introduced J to rotation rates measured by the G Ring laser in Germany. It is important to note that data are unavoidably affected by factors such as topography, instrument errors, background noise and so on. Furthermore, the bedrock at Wettzell consists mostly of igneous and metamorphic rocks, thus suggesting a very small correction. Nevertheless, we saw that the misfit was better between the observed and the corrected rotation rates as compared to the one between observed and synthetic rotation rates.

One other consequence of small-scale heterogeneities is that wavefield gradients measured at a point do not coincide with wavefield gradients estimated at the point using an array. However, we showed that the coupling tensor J for one can be related to coupling tensor J of another. This was confirmed in the last section where we see a disagreement between the rotation rates estimated at the G Ring laser using an array and the rotation rates measured by it. We corrected for the arrayderived rotation rates and the agreement between corrected rotation rates and the rotation rates measured by the ring laser improved significantly.

Although all examples here are based on rotation and rotation rates, the results can be extended to any other wavefield gradient with no modifications. Future work will focus on the relationship between J and subsurface geology. This relationship is highly non-linear and strong apriori information with Bayesian style inversion will be needed. Similarly, in section we perform the inversion for J for the ring laser at Wettzell, Germany. In this case, our effective model was PREM. Figure A2 shows that displacement data at the collocated seismometer WET and synthetics derived in this effective model are in agreement.

Therefore, the use of PREM here is justifiable. Note that the time axis are different than in Figure 7 to adjust the Rayleigh wave arrivals.

In section 4.2, we estimated array derived rotation rates at the ring laser. From Figure A3, we see that the velocities at three stations are in agreement. The distance between the stations is around 100 km.

APPENDIX D: SOURCE-RECEIVER DISTRIBUTION

Here, we show the map of the earthquakes and station data set used in 4. Figure A4 shows the source-receiver distribution for the inversion in 4.1, while Figure A5 shows the source-receiver distribution for the inversion in 4.2.

Figure 1 .

 1 Figure 1. Top: 2D domains used in section 2.1. Both domains are 20×20 km 2 . All boundaries are absorbing. The star shows the location of the source and the diamonds show the locations of four receivers. Left: Entirely homogeneous medium. Right: A heterogeneity of 200×200 m 2 , placed at the centre, is shown by the square (not to scale).

Figure 2 .

 2 Figure 2. Left: Vertical displacements (top), vertical velocities (second), vertical accelerations (middle), shear strains (fourth) and rotations (bottom) in the homogeneous medium. Right: Vertical displacements (top), vertical velocities (second), vertical accelerations (middle), shear strains (fourth) and rotations (bottom) in the medium with heterogeneity. The scaling factors relative to the displacements for each panel are shown on the top left and the two media are shown in Figure 1.

Figure 3 .

 3 Figure 3. Left: 2D domain used in section 3.1. The domain is 20×20 km 2 and the heterogeneity at the centre, not to scale, is 200×200 m 2 . The stars show the location of the sources and the diamond shows the location of the receiver. All boundaries are absorbing. Middle: Zoom in into the heterogeneity. Right:

Figure 4 .

 4 Figure 4. Rotations. Panels a, b, c: Three events used for inversion of the corrector. Panel d: The event is used as validation using this corrector. It is shown in blue to indicate that it has not been used for the inversion.

Figure 5 .

 5 Figure 5. Top: 3D domain used in section 3.2. The domain is 40×40×10 km 3 , and the heterogeneity at the centre, not to scale, is 1334×1334×222 m 3 . The stars show the location of the sources and the diamond shows the location of the receiver. All but top boundary are absorbing.

Figure 6 .

 6 Figure 6. Rotations. Panels a, b, c: Three events are used for inversion of the corrector. Panel d: The event is used as validation using this corrector. It is shown in blue to indicate that it has not been used for the inversion.

Figure 7 .

 7 Figure 7. Rotation Top: 201611131102A South Island, New Zealand. One of the 32 events that are used to invert for the coupling tensor J. Bottom: 201210280304A Queen Charlotte Islands. An event from the validation data set. It is shown in blue to indicate that this event has not been used for the inversion.

Figure 8 .

 8 Figure 8. Left: location of the receivers used for the finite difference. Right: Location of the seismometer ROTZ (top) and WET (bottom).

  A figure showing the sourcereceiver distribution in shown in Appendix D. Comparing the array-derived rotation rates to those observed at the ring laser, we find that the misfits are always above 0.75 even though the velocities from three seismometers are always in agreement [see Appendix C]. Therefore, we follow the

Figure 9 .

 9 Figure 9. Rotation rates. Top: 201103110546A Near East Coast of Honsh. One of the 6 events that are used to invert for the coupling tensor J. Bottom: 200909291748A Samoa Islands Region. An event from the validation data set. It is shown in blue to indicate that this event has not been used for the inversion.

Figure A1 .

 A1 Figure A1. Vertical displacements for the four events shown in Figure 4. Displacements in the heterogeneous media are shown in black and the displacements in the homogeneous media are shown in green.

Figure A2 .

 A2 Figure A2. Vertical displacements measured at WET for the two events shown in Figure 7. Observed displacements are shown in black and the displacements computed in PREM are shown in green. Top: 201611131102A South Island, New Zealand. Bottom: 201210280304A Queen Charlotte Islands.

Figure A3 .

 A3 Figure A3. Velocities along the East component for the two events shown in Figure 9. Observed velocities at stations ROTZ, WET and GRC1 are shown in black, green and red respectively. Top: 201103110546A Near East Coast of Honsh. Bottom: 200909291748A Samoa Islands Region.

Figure A4 .

 A4 Figure A4. The source receiver geometry for the events in Section 4.1. The 32 events used for the inversion are shown in red, while the event that was used for validation is shown in blue. The blue dot shows the location of the receiver.

Figure A5 .

 A5 Figure A5. The source receiver geometry for the events in Section 4.2. The 6 events used for the inversion are shown in red, while the event that was used for validation is shown in blue. The blue dots show the location of the three receivers.

ACKNOWLEDGEMENTS

We would like to thank Michael Afanasiev and Johana Brokešová whose comments have really helped improve the paper. This work was funded by the French National Research Agency (ANR-16-CE31-0022). We also thank the support through the European Research Council (Advanced Grant: ROMY, number: 339991).

Numerical simulations were performed at CCIPL and we thank the people who work endlessly to maintain the system. Obspy was acquired from https://github.com/obspy/obspy/wiki and all event data were acquired from the global CMT catalog (Ekström et al. 2012).

APPENDIX A: 2D NON-PERIODIC HOMOGENIZATION

Here, we give a brief overview on the theory of homogenization for the 2D case. For a more comprehensive description, Capdeville et al. (2010b), Capdeville et al. (2010a) and [START_REF] Guillot | 2-D non-periodic homogenization of the elastic wave Correcting wavefield gradients for the effects of local small-scale heterogeneities 27 equation: SH case[END_REF], [START_REF] Cupillard | Non-periodic homogenization of 3-D elastic media for the seismic wave equation[END_REF] are proposed for the 1D, 2D and 3D case respectively.

Homogenization deals with the obtaining homogenized or effective equations for systems with fine microscopic structure. The limit between the macroscopic and microscopic variables is denoted by ε = λ 0 λm , where λ 0 is the length-scale below which scales are considered microscopic and λ m is the minimum wavelength in the medium.

To separate the microscopic and macroscopic variables explicitly, a new space variable y = x ε is introduced such that x encapsulates all the macroscopic variations and y encapsulates all the microscopic variations. It is assumed that x and y can be treated as independent variables. Consequently, all partial derivatives with respect to x are redefined as

In homogenization, displacement u(x, t) and stress σ(x, t) are postulated to be two-scale asymptotic expansions in ε:

where the i are powers in ε but not on the u and σ. Substituting Equations (A.1) and (A.2) into the wave equation yields a series of equations which can be solved for each i. For i = 0, it emerges that u 0 does not depend on the microscopic variable y, i.e. to the leading order,

Unlike i = 0, the higher order terms in the series (A.2) do depend on y, for example,

where χ 1 is the first-order corrector, a third-order tensor. Going further, we also obtain the equation of the cell problem:

where c is the elastic tensor.

In fact, u 0 , the zeroth order effective displacement, is a solution of the classical wave equation

where the effective model (ρ * , c * ) can be obtained by solving the cell problem (Capdeville et al. 2010a).

The corresponding equations for gradients can be derived using (A.1), (A.3) and (A.4).

APPENDIX B: CORRECTOR FOR ARRAY-DERIVED WAVEFIELD GRADIENT

It is possible to estimate wavefield gradients just from displacements. The most straightforward procedure is based on truncated Taylor series,

where j ∈ {x 1 , x2 }. This approximation is correct to first-order in grid spacing h and is only valid if f is smooth with respect to h. Recalling the equation for true displacement from (7),

Expanding the terms in (B.2) individually, noting that the corrector χ is not smooth with respect to x, and using the relation y = x ε , we have,

2) and (B.3) in (B.1), we get,

u ε (x + hj, y + h ε j, t) -u ε (x, y, t) h = ∇ j u 0 (x, t) + J(y) : (u 0 (x, t))

APPENDIX C: TRANSLATIONS ARE IN GOOD AGREEMENT

As stated in section 2.3, the effective model needed for the inversion of J is not limited to the homogenized model. Any tomographic model that fits displacement data reasonably well can be used.

In section 3.1, we used the homogeneous medium as our effective model and it is indeed sensible to do so because as shown by Figure A1, the displacements in both the heterogeneous and homogeneous media are in agreement. This reasoning is also extendable to inversion in section 3.2.