Guan Huang
email: huangguan@tsinghua.edu.cn

Sergei Kuksin
email: sergei.kuksin@imj-prg.fr

ON AVERAGING AND MIXING FOR STOCHASTIC PDES

We examine the convergence in the Krylov-Bogolyubov averaging for nonlinear stochastic perturbations of linear PDEs with pure imaginary spectrum and show that if the involved effective equation is mixing, then the convergence is uniform in time.

Introduction

The Krylov-Bogolyubov averaging for stochastic PDEs which we are concerned with in this work, means the following. Starting with a linear PDE on a torus (or on a bounded domain) with pure imaginary spectrum we consider its ε-small nonlinear stochastic perturbation. Then the above mentioned averaging describes the behaviour of the distributions of actions of solutions for the perturbed equation on time-intervals of order ε -1 . Here the actions of solutions are made by the halves of squared norms of their Fourier coefficients with respect to the basis, made by eigenfunctions of the original linear system. Description of the limit is made via an auxiliary effective equation which is another nonlinear stochastic equation whose nonlinearity is made from resonant terms of the nonlinear part of the perturbation. The effective equation may be mixing, and then as time goes to infinity its solutions converge in distributions to a statistical equilibrium, given by a measure in a function space. An observation which we make in this work is that in the mixing case the convergence in distribution of actions of solutions for the perturbed equation to those of solution for the effective equation, described by the Krylov-Bogolyubov averaging, is uniform in time.

The Krylov-Bogolyubov averaging under discussion applies to various classes of stochastic PDEs, depending on the type of the original unperturbed linear system. In Sections 2-3 we discuss in details the averaging for stochastic complex Ginzburg-Landau (CGL) equations, regarded as perturbations of linear Schrödinger equations, and in Section 4 briefly repeat the argument for stochastic nonlinear wave equations.

Notation. For a Banach space B and R > 0 we denote BR (B) = {b ∈ B ∶ b B ≤ R}; for a metric space M , P(M ) stands for the space of probability Borel measures on M . By ⇀ we denote the weak convergence of measures and by D(ξ) we denote the distribution of a random variable ξ. For a function f and a measure µ we denote ⟨f, µ⟩ = ∫ f dµ.

CGL: the setting and result

We consider a stochastic CGL equation on a torus

T D ∶= R (L 1 Z) × R (L 2 Z) × ⋯ × R (L D Z), L 1 , . . . , L D > 0, u t + i(-∆ + V (x))u = εµ∆u + εP(∇u, u) + √ εη(t, x), u = u(t, x), x ∈ T D , (2.1)
where µ ∈ {0, 1}, P ∶ C D+1 → C is a C ∞ -smooth function, ε ∈ (0, 1] is a small parameter, the random force η(t, x) is white in time and regular in x, and the potential V (x) is a real smooth function. If µ = 0, the nonlinearity P(∇u, u) should not depend on ∇u. For simplicity we assume that µ = 1 (the case µ = 0 can be treated similarly). Again only to simplify presentation we also assume that V (x) > 0 for all x.

For any s ∈ R we denote by H s the Sobolev space of complex functions on T D , provided with the norm ⋅ s , u 2 s = ⟨(-∆) s u, u⟩ + ⟨u, u⟩, if s ⩾ 0, where ⟨⋅, ⋅⟩ is the real scalar product in L 2 (T D ; C),

⟨u, v⟩ = R T d uv dx, u, v ∈ L 2 (T D ; C).
Let {e l (x), l ∈ N} be the usual trigonometric basis of the space L 2 (T D ), parametrized by natural numbers. Then -∆e l = κ l e l , κ l ≥ 0, and we assume that 0 = κ 1 < κ 2 ≤ κ 3 , . . . . We take the force term η(t, x) in (2.1) to be of the form

η(t, x) = ∂ ∂t ξ(t, x), ξ(t, x) ∶= l⩾1 b l β l (t)e l (x). (2.2)
Here

β l (t) = β R l (t) + iβ I l (t)
, where β R l (t), β I l (t), l ⩾ 1, are independent real-valued standard Brownian motions, defined on a complete probability space (Ω, F, P) with a filtration {F t ; t ⩾ 0}. 1 As a function of x, ξ(t, x) is assumed to be smooth in the sense that the real numbers b l , l ≥ 1, decays to zero faster than any negative degree of l.

Introducing the slow time τ = εt, we rewrite eq. (2.1) as

u + ε -1 (-∆ + V (x))u = ∆u + P(∇u, u) + ξ(τ, x), u = u(τ, x), x ∈ T D , (2.3)
where the upper dot stands for ∂ ∂τ , and ξ(τ, x) is as in (2.2) with t ∶= τ and with another set of standard independent complex Brownian motions β l . Here and below we write stochastic PDEs with additive noise as nonlinear PDEs with forcing terms of the form (2.2).

Definition 2.1. If L and E are Banach spaces with norms ⋅ L and ⋅ E , then Lip m (L, E), m ≥ 0, is the collection of maps

F ∶ L → E such that for any R ⩾ 1, sup R>0 (1 + R ) -m Lip(F BR (L) ) + sup v∈ BR (L) F (v) E < ∞,
where Lip(f ) is the Lipschitz constant of a mapping f .

We make the following assumption concerning the well-posedness of eq. (2.3). There and everywhere below in our paper either always the indices s of involved Sobolev spaces H s are integer, or always they are any real numbers. Assumption 2.2. There exist numbers 0 < s 1 < s 2 < +∞ and m ∈ N such that for each s ∈ (s 1 , s 2 ),

(1) the mapping H s → H s-1 : u ↦ P(∇u, u) belongs to Lip m(H s , H s-1 );

(2) for any ε ∈ (0, 1] and u 0 ∈ H s equation (2.3) has a unique strong solution u ω (τ ; u 0 ), equal u 0 at τ = 0, defined for τ ⩾ 0, and

E sup θ⩽τ ⩽θ+1 u(τ ; u 0 ) 2 m s ⩽ C s ( u 0 s ), ∀θ ⩾ 0,
where C s (⋅) is a continuous and non-decreasing function.

(3) If we work in the category of Sobolev spaces with integer indices, then the integer segment (s 1 , s 2 ) ∩ Z contains at least two points.

1 So {β l (t)} are standard independent complex Brownian motions.

Under the above assumptions the family of solutions {u ω (τ ; u 0 ), u 0 ∈ H s } defines in the spaces H s , s ∈ (s 1 , s 2 ), Markov processes. Assumption 2.2 is satisfied for many nonlinearities P. In particular it holds if

P(∇u, u) = -u + zf p ( u 2 )u, z ∈ C, z = 1, Iz ≤ 0, Rz ≤ 0, (2.4)
where f p (r) is a non-decreasing smooth function on R, equal r p for r ≥ 1. The degree p ≥ 0 is any if D = 1, 2, and p < 2 (D -2) if D ≥ 3. See [4, Section 5]. We denote by A V the Schrödinger operator

A V u ∶= -∆u + V (x)u.
Let {λ l } l⩾1 be its eigenvalues, ordered in such a way that

0 < λ 1 ⩽ λ 2 ⩽ λ 3 ⩽ ⋯
(we recall that V > 0) and let {ϕ l , l ⩾ 1} ⊂ L 2 (T D ) be an orthonormal basis, formed by the corresponding eigenfunctions. We denote Λ = (λ 1 , λ 2 , . . . ) and call Λ the frequency vector of eq. (2.3). For a complex-valued function u ∈ H s we denote by

Ψ(u) ∶= v = (v 1 , v 2 , . . . ), v k ∈ C,
the vector of its Fourier coefficients with respect to the basis {ϕ l } l⩾1 : u = ∑ l⩾1 v l ϕ l . Note that Ψ is a real operator: it maps real functions u(⋅) to real vectors v. In the space of complex sequences v = (v 1 , v 2 , . . . ) we introduce the norms

v 2 s = ∞ k=1 ( λ k s + 1) v k 2 , s ∈ R ,
and denote h s = {v ∶ v s < ∞}. Then Ψ defines an isomorphism between the spaces H s and h s , for any s. Now we write Eq. (2.3) in the v-variables:

vk + ε -1 iλ k v k = -λ k v k + P k (v) + ∞ l=1 Ψ kl b l βl , k ∈ N. (2.5)
Here the k-th equation is obtained as the L 2 -scalar product with ϕ k of eq. (2.3), where u = Ψ -1 v. So

P (v) = (P k (v), k ∈ N) = Ψ V (x)u + P(∇u, u) , u = Ψ -1 v,
and Ψ kl = ⟨ϕ k , e l ⟩ (thus Ψ kl is the matrix of the operator Ψ with respect to the trigonometric basis in L 2 (T D ) and the natural basis in h 0 ). Our task is to study the dynamics of eq. (2.5) when ε ≪ 1. An efficient way to deal with this problem is through the interaction representation, which means transition from variables {v k (τ )} to variables {a k (τ )}, where

a k (τ ) = e iε -1 λ k τ v k (τ ), k ≥ 1.
In the a-variables equations (2.5) read

ȧk (τ ) = -λ k a k + e iε -1 λ k τ P k (Φ -ε -1 Λτ a) + e iε -1 λ k τ ∞ l=1 Ψ kl b l βl , k ∈ N , τ ≥ 0, (2.6)
where for a vector θ = (θ k , k ∈ N) ∈ R ∞ , Φ θ stands for the rotation in h s , defined by

Φ θ v = v ′ , v ′ k = e iθ k v k ∀ k . (2.7)
Clearly operators Φ θ define isometries of all spaces h s . By a ε (τ ; v 0 ) = a ε k (τ ; v 0 ), k ≥ 1) we denote a solution u(τ ; u 0 ), written in the a-variables. It solves system (2.6) with the initial data a(0) = v 0 ∶= Ψ(u 0 ), and in view of Assumption 2.2.( 2),

E sup θ⩽τ ⩽θ+1 a ε (τ ; v 0 ) 2 m s ⩽ C ′ s ( v 0 s ), ∀θ ⩾ 0. (2.8)
In order to describe the dynamics of eq. (2.6) with ε ≪ 1 we introduce an effective equation:

ȧk = -λ k a k + R k (a) + ∞ l=1 B kl βl , k ∈ N, τ ≥ 0.
(2.9)

Here {B kl , k, l ⩾ 1} is the principal square root of infinite matrix {A kl , k, l ⩾ 1},2

A kl = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∑ j⩾1 b 2 j Ψ kj Ψ lj , if λ k = λ l , 0, else,
and

R(a) ∶= (R k (a), k ∈ N) = lim T →∞ 1 T T 0 Φ Λt P (Φ -Λt a)dt (2.10)
(see [START_REF] Huang | Time-averaging for weakly nonlinear CGL equations with arbitrary potential[END_REF]Lemma 4] concerning this limit). Under (1) of Assumption 2.2, for s ∈ (s 1 , s 2 ), the mapping h s → h s-1 ∶ v ↦ R(v) belongs to Lip m(h s , h s-1 ), see [START_REF] Huang | Time-averaging for weakly nonlinear CGL equations with arbitrary potential[END_REF]. Therefore for any a 0 ∈ H s with s ∈ (s 1 , s 2 ), a strong solution a 0 (τ ; a 0 ) of eq. (2.9), equal a 0 at τ = 0, is unique and exists at least locally. In [4, Theorem 2 and Proposition 1] we proved the following result:

Theorem 2.3. If Assumption 2.2 holds, then for any s 1 < s * < s < s 2 and any v 0 ∈ h s we have: i) eq. (2.9) has a unique strong solution a 0 (τ ; v 0 ), τ ≥ 0, equal to v 0 at τ = 0. It belongs to C([0, ∞), h s) a.s., and for any θ ≥ 0,

E sup τ ∈[θ,θ+1] a 0 (τ ; v 0 ) 2 m s ⩽ C ′ s( v 0 s). (2.11)
ii) For any T > 0,

D a ε (τ ; v 0 ) [0,T ] ⇀ D a 0 (τ ; v 0 ) [0,T ] in P C([0, T ]); h s * ) as ε → 0. (2.12) Remark 2.4. 1)
In [START_REF] Huang | Time-averaging for weakly nonlinear CGL equations with arbitrary potential[END_REF] it was proved only the second assertion of the theorem, implying that a 0 is a solution of eq. (2.9) in h s * . Then (2.8) and Fatou's lemma imply (2.11). Indeed, for M ∈ N let Π M (a 1 , a 2 , . . . ) = (a 1 , . . . , a M , 0, . . . ). Then by (2.8) and convergence (2.12), for each θ ≥ 0 and any M, N ∈ N

E sup τ ∈[θ,θ+1] N ∧ Π M a 0 (τ ; v 0 ) 2 m s ⩽ C ′ s( v 0 s).
Sending first N → ∞ and then M → ∞ and using Fatou's lemma we recover (2.11). Since a 0 is a solution in h s * , then dues to (2.11) it is a solution in h s. Uniqueness of a solution is obvious.

2) It follows immediately that (2.12) also holds for a solution a ε (τ ; v 0 ) if the initial data v 0 is a r.v. in h s, independent from the random field ξ. Moreover, a simple analysis of the proof in [START_REF] Huang | Time-averaging for weakly nonlinear CGL equations with arbitrary potential[END_REF] implies that if v ω 0 s ≤ M a.s., then the rate of convergence (2.12) depends only on M (and, of course, on s * < s).

3) Theorem's assertion with the same proof remains true if in the r.h.s. of eq. ( 2.3) we replace the viscosity ∆u by the hyperviscosity -(-∆ + 1) r u, r ∈ N (provided that Assumption 2.2 holds for the equation). The equations with hyperviscosity are important for some applications, see below Example 2.11. 4) In [START_REF] Huang | Time-averaging for weakly nonlinear CGL equations with arbitrary potential[END_REF] the bounds E sup θ⩽τ ⩽θ+1 u(τ ; u 0 ) m s are assumed for all m-th moments of solutions, m ∈ N. It was done only for simplicity, and under the additional restriction in Assumption 2.2.(1), which specifies the growth of the the nonlinearity P, only the bounds for the moments with m ≤ 2 m are needed for the proof. This difference is rather insignificant since using the standard techniques of exponential supermartingales it is usually easy to obtain bounds for all m-th moments after the second moments are estimated.

Again, solutions a 0 (τ ; v 0 ) of (2.9) define in the space h s * a Markov process. Our goal in this work is to prove that if effective equation (2.9) is mixing, then convergence (2.12) is such that D(a ε (τ ; v 0 )) ⇀ D(a 0 (τ ; v 0 )) uniformly in τ ≥ 0, with respect to the dual-Lipschitz distance. We recall Definition 2.5. Let M be a complete and separable metric space. For any two measures µ 1 , µ 2 ∈ P(M ) the dual-Lipschitz distance between them is

µ 1 -µ 2 * L,M ∶= sup f ∈C(M ), f L,M ⩽1 ⟨f, µ 1 ⟩ -⟨f, µ 2 ⟩ ⩽ 2, where f L,M = Lipf + f C(M ) .
For future usage we note that for s ≥ s ′ space P(h s ) is naturally embedded in P(h s ′ ), and it easily follows from the definition that for µ 1 , µ 2 ∈ P(h s ) the distance

µ 1 -µ 2 * L,h s ′ is a non-decreasing function of s ′ ≤ s.
To proceed we have to assume that the effective equation is mixing, and the rate of mixing is uniform for initial data from bounded sets: Assumption 2.6. For some s * ∈ (s 1 , s 2 ) effective equation (2.9) is mixing in the space h s * with a stationary measure µ 0 ∈ P(h s * ), and for each M > 0 and v ∈ BM (h s * ), we have

D(a 0 (τ ; v)) -µ 0 * L,h s * ⩽ g M (τ ), (2.13)
where g is a continuous function of (M, τ ) which goes to zero when τ → ∞.

Relation (2.13) is a mild specification of the mixing in eq. (2.9), and a proof of the latter in fact usually establishes the former.

Our main result is the following:

Theorem 2.7. Under Assumptions 2.2 and 2.6, for any s ∈ (s * , s 2 ) and any v 0 ∈ h

s lim ε→0 sup τ ⩾0 D(a ε (τ ; v 0 )) -D(a 0 (τ ; v 0 )) * L,h s * = 0,
where a ε (τ ; v 0 ) and a 0 (τ ; v 0 ) solve respectively equations (2.6) and (2.9) with initial conditions a ε (0; v 0 ) = a 0 (0; v 0 ) = v 0 . Moreover, for any M > 0 the above convergence is uniform for v 0 ∈ BM (h s).

For v = (v 1 , v 2 , . . . ) ∈ h s we introduce the vector of action variables

I(v) = (I 1 (v), I 2 (v), . . . ), where I k (v) = 1 2 v k 2 , k = 1, 2, . . . . Then I(v) ∈ h s I ∩ R ∞ + , where h s I is the weighted l 1 -space with the norm I I,s = ∑ ∞ k=1 ( λ k s + 1) I k .
Since the interaction representation does not change the actions, then for the action variables of solutions for the original equations (2.5) we have Corollary 2.8. Under the assumptions of Theorem 2.7 actions of a solution v ε (τ ; v 0 ) for eq. (2.5)

with v 0 ∈ h s, s > s * satisfy lim ε→0 sup τ ⩾0 D I(v ε (τ ; v 0 )) -D (I(a 0 (τ ; v 0 )) * L,h s * I = 0. 2.1. Examples. For a finite-dimensional stochastic equation ȧk (τ ) = e iε -1 λ k τ P k (Φ -ε -1 Λτ a) + e iε -1 λ k τ N l=1 M kl βl (τ ), k = 1, . . . ,
⟨P(u) -P(v), u -v⟩ ≥ 0 ∀u, v ∈ L ∞ (T D ; C).
From definition (2.10) we see that operator R is monotone as well. Taking any two vectors v 1 , v 2 ∈ h s and denoting w(τ ) = a 0 (τ ; v 1 )a 0 (τ ; v 2 ) we immediately derive from the equation for w that

(d dt) w ω (τ ) 2 0 ≤ -1 2 λ 1 w ω (τ ) 2 0 . So w ω (τ ) 0 ≤ e -λ1τ v 1 -v 2 0 ∀ τ ≥ 0, ∀ω.
This very strong a-priori estimate implies Assumption 2.6 after some non-complicated work.

Example 2.10. Now let P(∇u, u) = P(u) has the form (2.4), where Rz < 0 and Iz < 0. Then relation (2.13) with g M (τ ) = C(M )e -κτ , κ > 0, and with sufficiently large integer s * follows from the abstract theorems in [10, Theorem 2.1] and [5, Theorem 3.1.7] since solutions a 0 inherit estimates for solutions a ε via convergence (2.12). We leave to the reader details of this derivation (which are not quite trivial).

Example 2.11. Consider eq. (2.3) with P(∇u, u) = -iρ u 2 u, with V = 0 and with a hyperviscosity instead of viscosity:

u -iε -1 ∆u = -(-∆ + 1) r u -iρ u 2 u + ξ(τ, x), x ∈ T D L = R D (LZ D ). (2.15)
Here r ∈ N, ρ > 0 is a scaling factor and ξ has the form (2.2), where all b l 's are nonzero. Now A V = -∆, so its eigenfunctions are exponents e i2πs⋅x , s ∈ L -1 Z D . The action variables I s (τ ) of solutions u(τ, x) are naturally parametrised by s ∈ L -1 Z D . The theory of wave turbulence (WT), among other things, examine the behaviour of the expectations of the actions EI s (τ ) under the limit of WT:

ε → 0, L → ∞,
when ρ is properly scaled with ε and L. See [START_REF] Nazarenko | Wave Turbulence[END_REF] and the introduction in [START_REF] Dymov | Formal expansions in stochastic model for wave turbulence 1: kinetic limit[END_REF]. For any dimension D solutions of eq. (2.15) with r sufficiently large in terms of D satisfy Assumption 2.2 by the same straightforward proof as for eq. ( 2.3) with P as in (2.4), and the equation also meets Assumption 2.6 with sufficiently large integer s * in view of the abstract theorems, mentioned in Example 2.10. So Corollary 2.8 applies to eq. (2.15) with ρ and L fixed, under the limit ε → 0 and imply that the expectations of the actions of solutions for (2.15) converge to those of solutions 3 In [2, Theorem 2.9] the uniform in time convergence as above is proved for a class of systems (2.14). The proof in [START_REF] Dymov | Nonequilibrium statistical mechanics of weakly stochastically perturbed system of oscillators[END_REF] is based on the observation that the corresponding v-equations are mixing with a rate of mixing, independent from ε.

for the corresponding effective equation, uniformly in time. The limit ε → 0 in eq. (2.15) and in similar equations is known in nonlinear physics as the limit of discrete turbulence. See [START_REF] Nazarenko | Wave Turbulence[END_REF], [START_REF] Kuksin | Resonant averaging for small-amplitude solutions of stochastic nonlinear Schrödinger equations[END_REF] and [3, Sections 1.2, 12.1].

3. Proof of Theorem 2.7

Below we always assume Assumptions 2.2 and 2.6. We fix s > s * and v 0 ∈ h s as in Theorem 2.7. We abbreviate ⋅ * L,h s * to ⋅ * L and a (τ ; v 0 ) to a (τ ). For any T ′ ≥ 0 we denote by a 0 T ′ (τ ) a weak solution of the effective equation (2.9) such that D(a 0 T ′ (0)) = D(a ε (T ′ )). Note that solution a 0 T ′ (τ ) depends on ε and that a 0 0 (τ ) = a 0 (τ ; v 0 ). The following lemma follows from Theorem 2.3, Assumption 2.2.(2) and Remark 2.4.2). Lemma 3.1. For any δ > 0 and T > 0 there exists

ε 1 = ε 1 (δ, T ) > 0 such that if ε ≤ ε 1 , then sup τ ∈[0,T ] D(a ε (T ′ + τ )) -D(a 0 T ′ (τ )) * L,h s * ⩽ δ 2, ∀ T ′ ⩾ 0. (3.1)
Since by Assumption 2.6, D(a 0 (τ ; 0)) → µ 0 in P(h s * ) as τ → ∞, then from estimate (2.11) and Fatou's lemma we derive that

⟨ v 2 m s , µ 0 ⟩ ⩽ C ′ s(0) ∶= C 0 . (3.2)
We need the following statement for the mixing in eq. (2.9).

Lemma 3.2. For any solution a 0 (τ ; µ) ∈ h s * , τ ⩾ 0, of effective equation (2.9) such that D(a 0 (0; µ)) =∶ µ, where the measure µ satisfies ⟨ a 2 m s * , µ(da)⟩ = E a 0 (0; µ) 2 m s * ⩽ M for some M > 0, we have

D a 0 (τ ; µ) -µ 0 * L ≤ g M (τ, d) if µ -µ 0 * L ≤ d ≤ 2. (3.3)
Here the function g ∶ R 3 + → R + , (M, τ, d) → g M (τ, d) is continuous, vanishes with d, converges to zero when τ → ∞ and is such that for each fixed M ⩾ 0 the function

(τ, d) → g M (τ, d) is uniformly continuous in d for (τ, d) ∈ [0, ∞) × [0, 2].
A proof of the lemma is rather straightforward but lengthy. It is given below in Subsection 3.1.

We denote

M * ∶= C ′ s( v 0 s), where C ′ s( v 0 s) is as in (2.8
) and (2.11). Constants in the estimates below depend on M * , but this dependence usually is not indicated. Lemma 3.3. Take any δ > 0 and choose a T * = T * (δ) > 0, satisfying

g M * (T, 2) ≤ δ 4 ∀ T ≥ T * .
Then there exists

ε 2 = ε 2 (δ) > 0 such that if ε ≤ ε 2 and D(a ε (T ′ )) -µ 0 * L ⩽ δ for some T ′ ≥ 0, then D(a ε (T ′ + T * )) -µ 0 * L ⩽ δ, (3.4)
and

sup τ ∈[T ′ ,T ′ +T * ] D(a ε (τ )) -µ 0 * L ⩽ δ 2 + sup τ ≥0 g M * (τ, δ). (3.5) Proof. Let us choose ε 2 (δ) = ε 1 ( δ 2 , T * (δ))
, where ε 1 (⋅) is as in Lemma 3.1. Then we get from (3.1), (2.11), (3.3) with M = M * and the definition of T * that for ε ≤ ε 2 ,

D(a ε (T ′ + T * )) -µ 0 * L ⩽ D(a ε (T ′ + T * )) -D(a 0 T ′ (T * )) * L + D(a 0 T ′ (T * )) -µ 0 * L ⩽ δ.
This proves (3.4). Next, in view of (3.1) and (3.3),

sup θ∈[0,T * ] D(a ε (T ′ + θ)) -µ 0 * L ⩽ sup θ∈[0,T * ] D(a ε (T ′ + θ)) -D(a 0 T ′ (θ)) * L + sup θ∈[0,T * ] D(a 0 T ′ (θ)) -µ 0 * L ⩽ δ 2 + max θ∈[0,T * ] g M * (θ, δ).
This implies (3.5).

We are now ready to prove Theorem 2.7.

Proof. [of Theorem 2.7.] Let us fix arbitrary δ > 0 and take some 0 < δ 1 ≤ δ 4. Below in the proof the functions ε 1 , ε 2 and T * are as in Lemmas 3.1 and 3.3. i) By the definition of T * = T * (δ 1 ), (2.11) and (3.3),

D a 0 T ′ (τ )) -µ 0 * L ≤ δ 1 4 ∀ τ ≥ T * , (3.6)
for any

T ′ ≥ 0. ii) By Lemma 3.1, if ε ≤ ε 1 = ε 1 ( δ1 2 , T * ) > 0, then sup 0≤τ ≤T * D a ε (τ )) -D a 0 (τ ; v 0 ) * L ≤ δ1 2 . (3.7)
In particular, in view of (3.6) with T ′ = 0,

D a ε (T * ) -µ 0 * L < δ 1 . (3.8)
iii) By (3.8) and (3.4) with δ = δ 1 and with T ′ = nT * , n = 1, 2, . . . we get inductively that

D a ε (nT * ) -µ 0 * L ≤ δ 1 ∀ n ∈ N, (3.9)
if ε ≤ ε 2 = ε 2 (δ 1 ). iv) Now by (3.9) and (3.5) with δ = δ 1 , for any n ∈ N and 0

≤ θ ≤ T * , D a ε (nT * + θ)) -µ 0 * L ≤ δ 1 2 + sup θ≥0 g M * (θ, δ 1 ), (3.10) if ε ≤ ε 2 (δ 1 ). v) Finally, if ε ≤ ε # (δ 1 ) = min ε 1 , ε 2 , then by (3.7) if τ ≤ T * and by (3.6)+(3.10) if τ ≥ T * we have that D a ε (τ )) -D a 0 (τ ; v 0 ) * L ≤ δ 1 + sup θ≥0 g M * (θ, δ 1 ) ∀ τ ≥ 0.
By Lemma 3.2, for M * fixed the function g M * (θ, d) is uniformly continuous in d and vanishes at d = 0. So there exists δ * = δ * (δ), which we may assume to be ≤ δ 4, such that if δ 1 = δ * , then g M * (θ, δ 1 ) ≤ δ 2 for every θ ⩾ 0. Then by the estimate above,

D a ε (τ )) -D a 0 (τ ; v 0 ) * L ≤ δ if ε ≤ ε * (δ) ∶= ε # (δ * δ) > 0,
for every positive δ. This proves the theorem's assertion.

3.1. Proof of Lemma 3.2. In this proof we write solutions of effective equation (2.9) as a(τ ) (rather than a 0 (τ )). i) At this step, for any non-random v 1 , v 2 ∈ BM (h s * ) we denote a ω j (τ ) = a ω (τ, v j ), j = 1, 2, and examine the distance D(a 1 (τ )) -D(a 2 (τ )) * L as a function of τ and v 1v 2 s * . We assume that v 1v 2 s * ≤ d for some d ≥ 0 and set w ω (τ ) = a ω 1 (τ )a ω 2 (τ ). We have ẇω = -A V w ω + R(a ω 1 ) -R(a ω 2 ). So by Duhamel's principle,

w ω (τ ) = e -A V τ w(0) + τ 0 e -A V (τ -t) (R(a ω 1 (t)) -R(a ω 2 (t)))dt.
Here w(0) s * ≤ d, and by (1) of Assumption 2.2

R(a ω 1 (τ )) -R(a ω 2 (τ )) s * -1 ≤ C w ω (τ ) s * X ω (τ ), X ω (τ ) = 1 + a ω 1 (τ ) m s * + a ω 2 (τ ) m s * . For θ > 0 the norm of the operator e -A V θ ∶ H s * -1 → H s * is bounded by χ(θ), where χ(θ) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ Cθ -1 2 , 0 < θ ≤ 1, Ce -cθ , θ > 1, with some C, c > 0. So if max 0≤t≤τ X ω (t) ≤ K, then w ω (τ ) s * ≤ de -λ1τ + KC τ 0 χ(τ -l) w ω (l) s * dl.
Applying Gronwall's lemma we derive from here that

w ω (τ ) s * ≤ d(1 + CKe c1K ), ∀ τ ≥ 0. (3.11)
for some positive constants C, c 1 . Denote Y (T ) = sup 0≤t≤T X ω (t) . By (2.11),

EY (T ) ⩽ 2(1 + C ′ s * (M ))(T + 1). For K > 0 let Ω K (T ) be the event {Y (T ) ≥ K}. Then P(Ω K (T )) ≤ 2(1 + C ′ s * (M ))(T +1)K -1
, and a 1 (τ )-a 2 (τ ) s * = w ω (τ ) s * satisfies (3.11) for ω ∉ Ω K (T ). From here we see that

if f ∈ C b (h s * ) is such that f C b (h s * ) ≤ 1 and Lip f ≤ 1, then E f (a 1 (τ )) -f (a 2 (τ )) ≤ 2P(Ω K (τ )) + d(1 + CKe c1K ) ⩽ 4(1 + C ′ s * (M ))(τ + 1)K -1 + d(1 + CKe c1K ) ∀ K > 0.
(3.12)

Let us denote by g 1 M (τ, d) the function in the r.h.s. above with K = ln ln d-1 ∨ 3 . This is a continuous function of (M, τ, d) ∈ R 3 + , vanishing when d = 0. Due to (2.13) and (3.12),

D(a(τ ; v 1 ))-D(a(τ ; v 2 )) * L = D(a 1 (τ )) -D(a 2 (τ )) * L ≤ (2g M (τ )) ∧ g 1 M (τ, d) ∧ 2 =∶ g 2 M (τ, d), (3.13) if v 1 , v 2 ∈ BM (h s * ) and v 1 -v 2 s * ≤ d, for any M, d > 0.
The function g 2 is continuous in the variables (M, τ, d), vanishes with d and goes to zero when τ → ∞ since g M (τ ) does.

ii) At this step we consider a solution a 0 (τ ; µ) =∶ a(τ ; µ) of (2.9) as in the lemma and examine the l.h.s. of (3.3) as a function of τ . For any K > 0 consider the conditional probabilities µ K = µ(⋅ BK (h s * )) and μK = µ(⋅ h s * ∖ BK (h s * )). Then

µ = A K µ K + ĀK μK , where A K = µ( BK (h s * )) and ĀK = P{ a(0) s * > K} ⩽ M K 2 m as E a(0) 2 m s * ≤ M . So D(a(τ, µ)) = A K D(a(τ ; µ K )) + ĀK D(a(τ ; μK )). (3.14)
In view of Assumption 2.6,

D(a(τ ; µ K )) -µ 0 * L = D(a(τ ; v)) µ K (dv) -µ 0 * L ≤ D(a(τ ; v)) -µ 0 * L µ K (dv) ≤ g K (τ ). Therefore due to (3.14), D(a(τ, µ)) -µ 0 * L ⩽ A K D(a(τ, µ K )) -µ 0 * L + ĀK D(a(τ, μK )) -µ 0 * L ≤ D(a(τ, µ K )) -µ 0 * L + 2 ĀK ≤ g K (τ ) + 2 M K 2 m
for K > 0 and τ ≥ 0. Let K 1 (τ ) > 0 be a continuous non-decreasing function such that K 1 (τ ) → ∞ and g K1(τ ) (τ ) → 0 as τ → ∞ (it exists since g K (τ ) is a continuous function of (K, τ ), going to 0 as τ → ∞ for each fixed K). Then choosing in the estimate above K = K 1 we get

D(a(τ ; µ)) -µ 0 * L ⩽ g K1(τ ) (τ ) + 2M K 1 (τ ) 2 m =∶ ĝM (τ ). (3.15)
Obviously ĝM (τ ) ≥ 0 is a continuous function on R 2 + , converging to 0 as τ → ∞. iii) Now we examine the l.h.s. of (3.3) as a function of τ and d. Recall that the Kantorovich distance between measures ν 1 , ν 2 on h s * is

ν 1 -ν 2 K = sup Lip f ≤1 ⟨f, ν 1 ⟩ -⟨f, ν 2 ⟩ ≤ ∞. Obviously ν 1 -ν 2 * L ≤ ν 1 -ν 2 K .
Since the 2 m-th moments of µ and µ 0 are bounded by M ∨ C 0 by (3.2) and the assumption on µ and since µ -

µ 0 * L ≤ d, then µ -µ 0 K ≤ C(M ∨ C 0 ) γ1 d γ2 ∶= d, γ 1 = 1 2 m , γ 2 = 2 m-1 2 m ; (3.16)
see [START_REF] Boritchev | One-dimensional turbulence and the stochastic Burgers equation[END_REF]Section 11.4] and [START_REF] Villani | Optimal Transport[END_REF]Chapter 7]. By the Kantorovich-Rubinstein theorem (see [START_REF] Villani | Optimal Transport[END_REF][START_REF] Boritchev | One-dimensional turbulence and the stochastic Burgers equation[END_REF]) there exist r.v.'s ξ and ξ 0 , defined on a new probability space (Ω ′ , F ′ , P ′ ), such that D(ξ) = µ, D(ξ 0 ) = µ 0 and

E ξ 1 -ξ 0 s * = µ -µ 0 K . (3.17)
Then using (3.13) and denoting by a st (τ ) a stationary solution of (2.9), D(a st (τ )) ≡ µ 0 , we have:

D(a(τ )) -µ 0 * L = Da(τ ; µ 0 ) -D(a st (τ )) * L ≤ E ω ′ D(a(τ ; ξ ω ′ )) -D(a(τ ; ξ ω ′ 0 )) * L ≤ E ω ′ g 2 M (τ, ξ ω ′ -ξ ω ′ 0 s * ), M = M ω ′ = ξ ω ′ s * ∨ ξ ω ′ 0 s * . As E ω ′ M 2 m ≤ 2(M ∨ C 0 ), then for any K > 0, P ω ′ (Q ′ K ) ≤ 2K -2 m(M ∨ C 0 ), Q ′ K = { M ≥ K} ⊂ Ω ′ . Since g 2 ≤ 2 and for ω ′ ∉ Q ′ K we have ξ ω ′ s * , ξ ω ′ 0 s * ≤ K, then D(a(τ )) -µ 0 * L ≤ 4K -2 m(M ∨ C 0 ) + E ω ′ g 2 K (τ, ξ ω ′ -ξ ω ′ 0 s * ) ∀ K > 0. For an r > 0 let us denote Ω ′ r = { ξ ω ′ -ξ ω ′ 0 s * ≥ r}.
Then by (3.17) and (3.16),

P ω ′ Ω ′ r ≤ dr -1 . So D(a(τ ))-µ 0 * L ≤ 4K -2 m(M ∨C 0 )+2 dr -1 +g 2 K (τ, r), for any K, r > 0. (3.18)
Let g 0 (l) be a positive continuous function on R + such that g 0 (l) → ∞ as l → +∞ in such a way that C ′ s * g 0 (l) (ln ln l) -1 2 ⩽ 2C ′ s * (0) for l ⩾ 3. With r = d1 2 and K = g 0 (r -1 ), we denote the r.h.s of (3.18) as g 3 M (τ, r) (so we substitute in (3.18) d = r 2 and K = g 0 (r -1 )). By (3.18) and the definition of g 2 (see (3.13)), we have

D(a(τ )) -µ 0 * L ≤ g 3 M (τ, r) ⩽ 4(g 0 (r -1 )) -2 m(M ∨ C 0 ) + 2r + 4(τ + 1) 1 + C ′ s * g 0 (r -1 )
ln ln(r -1 ∨ 3) -1 + rC ln ln(r -1 ∨ 3) exp c 1 ln ln(r -1 ∨ 3) .

As r → 0 the second and fourth terms converge to zero. By the choice of g 0 , the first term clearly converges to zero with r, so does the third term, which is ⩽ 8(1 + τ )(1 + C ′ s * (0))(ln ln(r -1 )) -1 2 for r ⩽ M is uniformly continuous in d, and the lemma is proved.

NLW: the setting and result

In this section we briefly discuss stochastic nonlinear wave (NLW) equations. Following [START_REF] Martirosyan | Exponential mixing for the white-forced damped nonlinear wave equation[END_REF][START_REF] Martirosyan | Local large deviations principle for occupation measures of the stochastic damped nonlinear wave equation[END_REF] we consider the following equations on a smooth bounded domain

D ⊂ R 3 : ∂ 2 t u + γ∂ t u -∆u = -γf (u) + γh(x) + √ γ η(t, x), x ∈ D, (4.1
) supplemented with the Dirichlet boundary condition on ∂D. Here γ ∈ (0, 1] is a small parameter, h(x) is a function in H 1 0 (D; R) and the nonlinearity f is C 2smooth. The random force η(t, x) is a white noise of time of the form

η(t, x) = ∂ ∂t ∞ j=1 b j β R j (t)e j (x).
Here {β R j (t), j ⩾ 1} is a sequence of independent standard real Brownian motions, {e j (x), j ⩾ 1} is an orthonormal basis in L 2 (D; R) composed of eigenfunctions of the Laplacian operator; that is -∆e j = λ j e j with 0 < λ 1 ⩽ λ 2 ⩽ . . . . The set {b j , j ⩾ 1} is a sequence of positive real numbers, satisfying

B ∶= j≥1 λ j b 2 j < +∞. (4.2)
The nonlinear term f meets the following growth condition,

f ′′ (u) ⩽ C( u ρ-1 + 1), u ∈ R, (4.3)
where C and ρ < 2 are positive constants, as well as the dissipativity conditions

F (u) ⩾ -κu 2 -C, u ∈ R, (4.4)
f (u)u -F (u) ⩾ -κu 2 -C, u ∈ R, (4.5
) where F is the primitive of f and κ is a positive constant.

f (u) -f (v) 2 0 ⩽ C( u 2ρ 1-sρ + v 2ρ 1-sρ + 1) u -v 2 1-sρ , s ρ = 2-ρ 2(ρ+1) .
Let us denote ∂ t u = Lv, where L = (-∆) 1 2 , and set ξ = u + iv. Then in terms of ξ equation (4.1) reads

∂ t ξ = iLξ -γiIξ -γiL -1 f (Rξ) -γiL -1 η(t, x).
Introducing the slow time τ = γt, we have

∂ τ ξ = γ -1 iLξ -iIξ -iL -1 f (Rξ) -iη(τ, x), (4.6)
where η(τ, x) = ∑ ∞ j=1 bj βR j (τ )e j (x) with bj = b j λ -1 2 j

, j ⩾ 1. The NLW equation (4.1) (without the γ-factor in the r.h.s) was studied in [START_REF] Martirosyan | Exponential mixing for the white-forced damped nonlinear wave equation[END_REF][START_REF] Martirosyan | Local large deviations principle for occupation measures of the stochastic damped nonlinear wave equation[END_REF], where the global well-posedness and estimates for the norms of solutions are established in the Sobolev space

H s ×H s-1 = {(u, u)}, s ∈ [1, 2-ρ 2 ]
. A simple analysis of the proof of [START_REF] Martirosyan | Exponential mixing for the white-forced damped nonlinear wave equation[END_REF]Proposition 3.4] or [START_REF] Martirosyan | Local large deviations principle for occupation measures of the stochastic damped nonlinear wave equation[END_REF]Proposition 5.4] implies the following statement on the global well-posedness of eq. (4.6): Theorem 4.2. Assume conditions (4.3), (4.4) and (4.5). Then for any s ∈ [1, 2-ρ 2 ) and ξ 0 ∈ H s (D; C), equation (4.6) has a unique strong solution ξ γ (τ ; ξ 0 ), equal to ξ 0 at τ = 0 and defined for τ ⩾ 0. Uniformly in γ ∈ (0, 1] this solution satisfies

E sup τ ∈[θ,θ+1] ξ γ (τ ; ξ 0 ) 2m s ⩽ C s (m, ξ 0 s , B), ∀θ ⩾ 0, m ∈ N, (4.7)
where C s is a continuous function, non-decreasing in all its arguments. Now let us write eq. (4.6) in terms of Fourier coefficients with respect to the basis {e j (x)}: ξ = ∑ j⩾1 v j e j . As in Section 2, we use h s to denote the space of sequences of complex Fourier coefficients and denote by Ψ the map

H s (D; C) → h s , ξ ↦ v ∶= (v j , j ⩾ 1). Then we have vk = γ -1 iλ 1 2 k v k -iIv k + i Pk (v) + i bk β R k (τ ), k ∈ N, (4.8)
where

P (v) = ( Pk (v), k ∈ N) = Ψ(L -1 f (Rξ)) with ξ = Ψ -1 v. Passing to the inter- action representation a k (τ ) = e -iγ -1 λ 1 2 k τ v k (τ ), k ⩾ 1
, we obtain the following system of equations for the a-variables

ȧk (τ ) = ie -iγ -1 λ 1 2 k τ -I e iγ -1 λ 1 2 k τ a k + Pk (Φ γ -1 Λτ a) + bk βR k , k ∈ N, (4.9)
where Λ = (λ 1 2

k , k ⩾ 1) and operator Φ is defined in (2.7). Let us calculate an effective equation for (4.9), following [START_REF] Huang | Time-averaging for weakly nonlinear CGL equations with arbitrary potential[END_REF].

i) To calculate the effective diffusion we decomplexify components a k = a R k + ia I k of the a-vector as (a R k + a I k ) ∈ R 2 and write the dispersion matrix of the equation as a block-diagonal real matrix with the blocks bk cos ϕ ksin ϕ k 0 0 , where

ϕ k = ϕ k (τ ) = γ -1 λ 1 2
k τ for k = 1, 2, . . . . So the diffusion matrix of eq. (4.9) -let us call it A(τ ) -is formed by blocks b2

k cos 2 ϕ k -cos ϕ k sin ϕ k -cos ϕ k sin ϕ k sin 2 ϕ k . Consider the limit A(τ ) = lim γ→0 1 τ ∫ τ 0 A(l) dl
. This is a τ -independent block-diagonal real matrix with the blocks 1 2 b2

k id, and this is the diffusion matrix of the effective equation. Corresponding dispersion matrix has blocks 1 √ 2 bk id. Coming back to the complex coordinates we see that the noise in the effective equation is 1 √ 2 ∑ bk βk (τ ), where { βk (τ )} are standard independent complex Brownian motions.

ii) The drift in the effective equation is a sum of two terms. The second one is In particular, the mapping a ↦ R(a) belongs to Lip ρ+1 (h s , h s ), s ∈ [s ρ , 1]. Therefore effective equation (4.10) is at least locally well-posed in h s , s ∈ [s ρ , 1]. Let us denote a 0 (τ ; a 0 ) its solution with an initial condition a 0 (0; a 0 ) = a 0 . Similarly let a γ (τ ; v 0 ) be a solution of eq. (4.9), equal a 0 at τ = 0. Then using (4.7) and arguing as in [START_REF] Huang | Time-averaging for weakly nonlinear CGL equations with arbitrary potential[END_REF] we get where C is a continuous function, increasing in all its arguments. ii) For any s > 1, a 0 ∈ h s and T > 0 we have D a γ (τ ; a 0 ) [0,T ] ⇀ D a 0 (τ ; a 0 ) [0,T ] in P C([0, T ]); h 1 ) as γ → 0. (4.13)

R(a) = ( Rk (a), k ⩾ 1) ∶= lim T →∞ 1 T T 0 Φ -Λt P (Φ Λt a)dt (cf. (2.10)). The k-th component of the first term is lim T →∞ 1 T T 0 e -iλ 1 2 k t -iI(e iλ 1 2 k t a k ) dt = -i 1 2 (Ia k -iRa k ) = -
Moreover, for any M > 0 the above convergence is uniform for a 0 ∈ BM (h s ).

Repeating the argument in Section 3 we get the following analogy to Theorem 2.7: Moreover, for any M > 0 the above convergence is uniform for v 0 ∈ BM (h s ).

For the original equation (4.1) we introduce the vector of action variables for its solution u = (u, u t ) as I(u) = (I 1 (u), . . . ), where I k (u) = 1 2 ⟨u, e k ⟩ 2 + 1 2 λ -1 k ⟨u t , e k ⟩ 2 , k = 1, . . . . Then we have the following corollary in analogy to Corollary 2.8: The result of Theorem 4.4 is conditional since it requires that eq. (4.9) is mixing. It is not our goal in this short section to check the mixing property. But we mention that since the effective equation (4.10) is similar to the NLW equation (4.1), written in the form (4.6), and since by (4.13) solutions for (4.10) inherit the estimates on solutions for (4.1), obtained in [START_REF] Martirosyan | Exponential mixing for the white-forced damped nonlinear wave equation[END_REF][START_REF] Martirosyan | Local large deviations principle for occupation measures of the stochastic damped nonlinear wave equation[END_REF], then most likely the proof of the exponential mixing in eq. (4.1), given in [START_REF] Martirosyan | Exponential mixing for the white-forced damped nonlinear wave equation[END_REF] (and based on an abstract theorem from [START_REF] Kuksin | Mathematics of Two-Dimensional Turbulence[END_REF]), applies to establish the exponential mixing for eq. (4.10) and thus verify for it the analogy of Assumption 2.6, required in Theorem 4.4.

1 3 .

 3 Hence g 3 M (τ, r) defines a continuous function on R 3 + , vanishing with r. Using (3.16) let us write r = d1 2 as r = R M (d), where R is a continuous function R 2 + → R + , non-decreasing in d and vanishing with d. Setting g 4 M (τ, d) = g 3 M (τ, R M (d ∧ 2)) and using that µµ 0 * L ≤ 2, we get from the above that D(a(τ ))µ 0 * L ≤ g 4 M (τ, µµ 0 * L ). Finally, evoking (3.15) we arrive at (3.3) with g M = g 5 M , where g 5 M (τ, d) = g 4 M (τ, d) ∧ ĝM (τ ) ∧ 2, 0 ≤ d ≤ 2. The function g 5 is continuous, vanishes with d and converges to zero as τ → ∞. For any fixed M > 0 this convergence is uniform in d due to the term ĝM (τ ). So for a fixed M > 0 the function (τ, d) ↦ g 5 M (τ, d) extends to a continuous function on the compact set [0, ∞] × [0, 2] (where it vanishes when τ = ∞). Thus g 5

Lemma 4 . 1 . ([ 7 ,

 417 Lemma 4.5]). Under the condition (4.3) we have

1 2 1 √ 2 2 ) 1 ⩽ C a 1 ρ 1 -sρ ∧ a 2 ρ 1 -

 112211121 a k . iii) So the effective equation for (4.9) reads ȧk (τ ) = -1 2 a k + i Rk (a) + bk d βk , k ∈ N. (4.10) By Lemma 4.1 and [4, Lemma 4], we have R(a 1 ) -R(a sρ + 1 a 1a 2 1-sρ , a 1 , a 2 ∈ h 1 . (4.11)

Theorem 4 . 3 .

 43 Under the assumptions of Theorem 4.2 we have: i) for any a 0 ∈ h 1 , eq. (4.10) has a unique strong solution a 0 (τ ; a 0 ), τ ⩾ 0. It belongs to C([0, ∞), h 1 ) a.s., and for any θ ⩾ 0 satisfiessup τ ∈[θ,θ+1] a 0 (τ ; a 0 ) 2m 1 ⩽ C(m, a 1 , B), ∀m ∈ N,(4.12)

Theorem 4 . 4 .

 44 Assume in addition to the assumption in Theorem 4.3 that the effective equation (4.10) is mixing in the space h 1 and that the requirement as in Assumption 2.6 is met. Then for any s > 1 and v 0 ∈ h s we have the following convergence lim γ→0 sup τ ⩾0 D(a γ (τ ; v 0 )) -D(a 0 (τ ; v 0 )) * L,h 1 = 0.

Corollary 4 . 5 .L,h 1 I

 451 Under the assumption of Theorem 4.4, for any s > 1 and u 0 = [u 1 , u 2 ] ∈ H s × H s-1 , the action-vector of a solution u γ (t; u 0 ) for equation (4.1), equal u 0 at t = 0, satisfieslim γ→0 sup τ ⩾0 D I(u γ (τ γ -1 ; u 0 )) -D (I(a 0 (τ ; v 0 )) * = 0,where v 0 = Ψ(u 1 + iL -1 u 2 ).

  , . . . , a N )(τ ) ∈ C N and Λ = (λ 1 , . . . , λ N ), satisfying a natural analogy of Assumptions 2.2, 2.6, a natural version of Theorem 2.7 holds. If the matrix (M kl ) is non-degenerate and Assumption 2.2 is met, then a convenient sufficient condition for Assumption 2.6 follows from the Khasminski criterion for mixing (details will be given in paper on the finite-dimensional stochastic averaging under preparation).3 For infinite-dimensional stochastic systems an instrumental criterion of mixing is not known yet. We believe that as in finite dimensions, Assumption 2.6 holds for effective equations for various equations (2.3), satisfying Assumption 2.2, at least if the random forces (2.2) is such that all b l 's are non-zero. Examples below are given to support this belief.

	N,	(2.14)
	where a(τ ) = (a 1

Example 2.9. Consider eq. (2.3), where P(∇u, u) = P(u) has the form (2.4) with z = -1. Since the function f p is monotone, then mapping P also is monotone in the sense that

The matrix (A kl ) defines a non-negative compact self-adjoint operator in the space l 2 . So its principal square root (which defines another non-negative compact self-adjoint operator) exists.

Acknowledgment

We are thankful to Armen Shirikyan for discussion. GH was supported by National Natural Science Foundation of China (Significant project No.11790273).


Junliang Dong

Xiaolong Wu

Alexandre Locquet

David S Citrin

Terahertz Superresolution Stratigraphic Characterization of Multilayered Structures Using Sparse Deconvolution

published or not. The documents may come   L'archive ouverte pluridisciplinaire


Junliang Dong

Alexandre Locquet

D S Citrin

Depth-resolution enhancement of terahertz deconvolution by autoregressive spectral extrapolation

This Letter presents a method for enhancing the depth-resolution of terahertz deconvolution based on autoregressive spectral extrapolation. The terahertz frequency components with high signal-to-noise ratio are modeled with an autoregressive process, and the missing frequency components in the regions with low signal-to-noise ratio are extrapolated based on the autoregressive model. In this way, the entire terahertz frequency spectrum of the impulse response function, corresponding to the material structure, is recovered. This method, which is verified numerically and experimentally, is able to provide a 'quasi-ideal' impulse response function, and therefore, greatly enhances the depth-resolution for characterizing optically thin layers in the terahertz regime.

Due to the penetrative capability of terahertz (THz) radiation, THz imaging has attracted considerable interest as a noninvasive, noncontact, and nonionizing modality to characterize various non-metallic multi-layered structures in industrial [START_REF] Su | Terahertz sensor for non-contact thickness and quality measurement of automobile paints of varying complexity[END_REF][START_REF] Dong | Nondestructive evaluation of forced delamination in glass fiber-reinforced composites by terahertz and ultrasonic waves[END_REF][START_REF] Lin | Quantifying pharmaceutical film coating with optical coherence tomography and terahertz pulsed imaging: an evaluation[END_REF], biomedical [START_REF] Parrot | Terahertz pulsed imaging in vivo: measurements and processing methods[END_REF] and cultural-heritage fields [START_REF] Schwerdtfeger | Beating the wavelength limit: three-dimensional imaging of buried subwavelength fractures in sculpture and construction materials by terahertz time-domain reflection spectroscopy[END_REF][START_REF] Koch-Dandolo | Reflection terahertz time-domain imaging for analysis of an 18th century neoclassical easel painting[END_REF][START_REF] Dong | Terahertz frequency-wavelet domain deconvolution for stratigraphic and subsurface investigation of art painting[END_REF]. THz reflective imaging provides information in depth by analyzing the reflected THz signal with an incident short THz pulse. Due to dielectric discontinuities with depth, reflected temporal THz echoes associated with Fresnel coefficients between various interfaces are recorded as a function of transverse position in amplitude and time delay. The stratigraphy of multi-layered structures can be obtained based on the precise estimation of the arrival times of echoes in the reflected signal. However, when the layers of the structures are optically thin in the THz regime, the echoes may partially or totally overlap in time; therefore, when no prior knowledge of the structures is available, THz deconvolution is essential to resolve the overlapping echoes and extract the arrival times of superimposed echoes from noisy measurements.

Theoretically, the reflected THz signal r(t) is the convolution of the incident THz pulse i(t) with the impulse response function h(t), which should ideally consist of M time-shifted ideal impulses, with the number M corresponding to the material structure; therefore, the Fourier transform of this ideal impulse response function is broadband and should be a sum of M complex sinusoids throughout the entire frequency spectrum. Conventional deconvolution based on direct inverse filtering aims at retrieving the impulse response function by applying the inverse Fourier transform of the transfer function H inv ( f ), which is the ratio of the reflected to the incident THz spectra,

h inv (t) = FFT -1 H inv ( f ) = FFT -1 FFT r(t) FFT i(t) , (1)
where FFT denotes the Fourier transform and FFT -1 the inverse Fourier transform, and h inv (t) is the impulse response function obtained by directly inverse filtering. Frequently, successful deconvolution cannot be expected by directly applying Eq.1. Because the THz pulse is band-limited, with a fast-decreasing spectrum outside of the passband (in our case extending from 60 GHz to 3 THz), the division of small values due to the deficiency of THz sources in the low-and high-frequency regions, will give rise to large abnormal values as high spikes and ruin the signal-to-noise ratio there; therefore, H inv ( f ), as obtained from Eq.1 is only valid within the bandwidth with high SNR. Furthermore, these high spikes in the low-and high-frequency regions, where the SNR is low, will introduce severe low-and high-frequency noises in the impulse response function. THz frequency-wavelet domain deconvolution (FWDD) [START_REF] Chen | Frequency-wavelet domain deconvolution for terahertz reflection imaging and spectroscopy[END_REF] is specifically designed to enhance the deconvolution process by first employing frequency-domain filtering and then further improving the SNR by wavelet de-noising.

Baseline subtraction is also required to cancel the slow fluctuations corresponding to the low-frequency noise. Wiener filtering [START_REF] Chen | Frequency-wavelet domain deconvolution for terahertz reflection imaging and spectroscopy[END_REF], double Gaussian filter [START_REF] Chen | Total variation deconvolution for terahertz pulsed imaging[END_REF] and van Hann filter [START_REF] Dong | Terahertz frequency-wavelet domain deconvolution for stratigraphic and subsurface investigation of art painting[END_REF][START_REF] Galvão | Optimization of apodization functions in terahertz transient spectrometry[END_REF] can serve as the frequency-domain filtering. However, the frequency-domain filtering will narrow the bandwidth of the impulse response function by eliminating the high spikes in the low-and high-frequency regions, and only the frequency components within the bandwidth with high SNR are kept. Therefore, pulses in the impulse response function recovered by FWDD are much wider than the ideal impulses, and consequently, limit the depth-resolution.

In this Letter, we demonstrate a method for enhancing the depth-resolution of THz deconvolution based on autoregressive spectral extrapolation. This method consists in modeling part of the deconvolved spectrum with high SNR with an autoregressive (AR) process, and to extrapolate the remaining part of the spectrum based on this AR model. What underlies AR spectral extrapolation is the maximum entropy method which consists in determining a spectral estimation that maximizes uncertainty with respect to the unknown information, thus eliminating unreasonable constraints, but that is consistent with the known information [START_REF] Ulrych | Time series modelling and maximum entropy[END_REF]. Direct use of the maximum entropy method is marred, however, by the lack of criterion for determining the order of the model used. It has been shown [START_REF] Ulrych | Maximum entropy spectral analysis and autoregressive decomposition[END_REF] that the maximum entropy method is equivalent to the least squares fitting of an AR process for which mathematical criterion (e.g. Akaike's criterion [START_REF] Akaike | A new look at the statistical model identification[END_REF]) exist to determine the length of an AR model. The existence of such a criterion is crucial to avoid over-or under-fitting. In the following, we focus on deconvolution exploiting AR spectral extrapolation, and call this method 'AR deconvolution' for convenience. Unlike FWDD, which discards the frequency components in the low-SNR regions in order to cancel the high spikes, AR deconvolution aims at finding the missing frequency components in the low-SNR regions using an AR model based on the data in the high-SNR regions. With AR deconvolution, the spectrum of the impulse response function throughout the entire frequency band can be estimated, then the inverse Fourier transform leads to a 'quasi-ideal' impulse response function can be achieved, and therefore, the resolution in depth can be enhanced.

In an AR process, each data point of H i of a digitized signal is a weighted sum of its previous p points, plus a noise term N i ,

H i = p ∑ k=1 a k H i-k + N i , i > p, (2)
where p is the order of the AR process, a k are the AR coefficients, and H i is known within the frequency band (i L ⩽ i ⩽ i H ) with high SNR. This AR model can be used as a prediction filter to find an estimate Ĥi of the unknown values of H i for i > i H using the forward prediction equation,

Ĥi = p ∑ k=1 a k H i-k , i > i H . (3)
Similarly, we can use a backward prediction filter to find the values of

H i for i < i L , Ĥi = p ∑ k=1 b k H i+k , i < i L , (4)
where b k are the coefficients for the backward prediction filter. The optimum values for the coefficients a k and b k are determined by minimizing the squared error between the model and the N BW (N BW = i Hi L + 1) available data points with high SNR.

As mentioned above, the frequency spectrum of the ideal impulse response function, of a multi-layered system, should be the superposition of a limited number of complex sinusoids, each corresponding to a layer interface. It has already been shown that [START_REF] Ulrych | Time series modelling and maximum entropy[END_REF], for a noiseless case, this kind of signal, which consists of a sum of complex sinusoids, can be modeled as an AR process with an order equal to twice the number of sinusoids; for a noisy case, this kind of signal should be modeled as an AR process with an order much higher than the number of sinusoids. Akaike's Information Criterion (AIC) [START_REF] Akaike | A new look at the statistical model identification[END_REF], which is based on the principle of entropy maximization in information theory, provides a measure for the selection of the model order p, and the AIC of order p can be calculated using the residual sums of squares from regression S 2 p , (AIC

) p = N ln(S 2 p ) + 2p, ( 5
)
where N is the number of data points, and ln is the natural logarithm. AIC deals with the trade-off between the goodness of fit of the model and the complexity of the model, and according to Akaike's theory, the most accurate model has the smallest AIC value. For the THz frequency spectrum, both the forward and backward prediction filters are needed to find the missing data in the low-and high-frequency regions, and the AR coefficients are obtained by minimizing the sum of the forward and backward prediction errors ϵ 2 , (3) data within [i PH , i H ] are used for the fitting to determine the a k . The Burg method is utilized to minimize this error term [START_REF] Shakibi | Resolution enhancement of ultrasonic defect signals for crack sizing[END_REF]. By adding the Levinson-Durbin constraint [START_REF] Ortigueira | Global versus local minimization in least-squares AR spectral estimation[END_REF], the Burg method enables the AR coefficients to be determined by a fast recursive algorithm, and guarantees a stable prediction filter, which is important to estimate the spectrum throughout the entire frequency band.

ϵ 2 = i PL ∑ i=i L |H i - p ∑ k=1 b k H i+k | 2 + i H ∑ i=i PH |H i - p ∑ k=1 a k H i-k | 2 , (6)
A standard THz time-domain spectroscopy (TDS) system (Teraview TPS Spectra 3000) is employed in this study to perform THz reflective imaging at almost normal incidence in order to obtain the data used in the following numerical and experimental study. The THz reference signal f 0 [n], with n as the index of data points, is recorded by setting a metal plate at the sample position, shown in the inset of Fig. 1. The data sampling period in the measurement is set to T s = 0.0116 ps. Each recorded reflected temporal THz waveform contains 4096 data points, and the signal is averaged over 10 shots. With this setting, the entire frequency spectrum obtained by Fourier transform is from 0 THz to 85.99 THz, which also contains 4096 data points.

Numerical simulations are first performed to verify the AR deconvolution. An ideal impulse response function h 0 [n], which contains 4096 data points with the same sampling period T s , is assumed and shown in Fig. 1. h 0 [n] represents a simple two-layered structure, and the time interval between the peaks are 25T s and 125T s , corresponding to the thickness of each layer. h 0 [n] is convolved with the reference signal f 0 [n] to simulate the received signal r 0 [n]. e[n] represents the additive white noise, with which the SNR of r 0 is set to be 10 dB in this simulation.

The assumed frequency spectrum H 0 , which ideally consists of a sum of three complex sinusoids and is not bandlimited throughout the frequency band (from 0 THz to 85.99 THz), is obtained by the Fourier transform of h 0 . In Fig. 2(a), H 0 is compared with H inv , which is the direct division of the spectra of r 0 and f 0 . By comparison, the low-SNR frequency bands of H inv which contain high spikes corresponding to low-and high-frequency noises, and the high-SNR frequency band can be well separated and clearly identified. The THz frequency spectrum has a large bandwidth with high SNR, which enables us to establish an AR model with a high order. Considering the THz bandwidth and the SNR of our TDS system, we fit the AR model to the N BW available data with i L = 0.17 THz and i H = 1.5 THz. By setting i PL = 0.3 THz and i PH = 1.2 THz, we allocate the data within [0.3, 1.5] for the forward prediction filter, and the data within [0.17, 1.2] for the backward prediction filter. AIC values with different orders are calculated and shown in the inset of Fig. 3. The smallest AIC value is achieved when p equals to 44. Fig. 2(b) shows the detailed data allocation for this AR process. After fitting this AR model to the N BW available data, we apply this AR model as a prediction filter to extrapolate the data to find the missing data in the regions below 0.18 THz and above 1.5 THz (up to 85.99 THz). Note that this is well outside the bandwidth of our imaging system; the ability to recover data out to this high frequency depends on our having sufficient knowledge of the frequency spectrum to construct an accurate AR model. The estimated frequency spectrum H AR based on this AR model is shown in red in Fig. 2. By comparison between H AR and H 0 , we can find out that the estimated frequency components in H AR match the assumed data in H 0 quite well, which verifies the effectiveness of this AR model to recover the missing data. By simply performing the inverse Fourier transform of H AR , the deconvolved signal h AR , can be achieved.

In Fig. 3, h AR is compared with h FWDD , which is the deconvolved signal obtained by FWDD, and the assumed impulse response function h 0 . For FWDD, van Hann filter is selected to serve as the frequency-domain filter to suppress the high-frequency noise, and the maximum value of the cutoff frequency we can set is f c = 3.5 THz in order to obtain a satisfactory SNR of the deconvolved signal. More details about the FWDD algorithm we use can be found in [START_REF] Dong | Terahertz quantitative nondestructive evaluation of failure modes in polymer-coated steel[END_REF]. With FWDD, only the first and the third interface can be identified, and the first assumed layer with time interval 25T s (corresponding to an air gap with thickness of 43.5 µm) cannot be resolved. In contrary, h AR is able to resolve all the interfaces, and exhibits a 'quasi-ideal' impulse response function compared with the assumed impulse response function h 0 . Based on this simulation, we can conclude that the minimum thickness resolved by THz reflective imaging, is enhanced by AR deconvolution, since the frequency components throughout the entire spectrum are estimated and recovered. It is important to note that the minimum thickness (time interval) resolution , corresponding to the depth-resolution, is mainly dependent on the bandwidth of the THz source, the sampling frequency and the SNR. In our simulations, the minimum time interval which can be resolved by AR deconvolution with the SNR = 10 dB is 15T s (corresponding to an air gap with thickness of 26.1 µm); on the other hand, the minimum time interval resolved by FWDD with the same SNR is 30T s . Therefore, the depth-resolution can be increased by a factor of two using AR deconvolution. Next, a one-layered polymer coating on a steel substrate is employed in this study to verify the effectiveness of AR deconvolution experimentally. The thickness of the polymer coating is about 22.5 µm, which is thin in the THz regime. THz reflective imaging experiment was performed, and the typical waveform received is shown in the inset of Fig. 4(a). Since the coating is very thin and the steel substrate provides total reflection, the reflected signal is similar to the reference signal and no overlapping echoes can be observed in time domain. Both FWDD and AR deconvolution are performed to resolve this thin coating layer. In the deconvolution, we consider the THz reference signal as the input and the received THz signal as the output; therefore, the actual impulse response function associated with the Fresnel coefficients can be obtained by multiplying the deconvolved signal by a factor of -1 for phase correction. For FWDD, the maximum cutoff frequency leading to a satisfactory SNR in time domain is 4 THz; For AR deconvolution, we set the parameters the same as the ones we used in the simulation. Based on the AR model, the frequency components outside the frequency In this Letter, AR deconvolution which is based on autoregressive spectral extrapolation, is presented to enhance the depth-resolution of THz deconvolution. The de-noising procedure in the conventional deconvolution, such as FWDD, narrows the bandwidth of the impulse response function by eliminating the frequency components in the low-SNR regions, and therefore, lowers the resolution in-depth. On the contrary, AR deconvolution aims at recovering the missing frequency components in the low-SNR regions using the AR model based on the frequency components in the high-SNR regions. In this way, the entire spectrum of the impulse response function is estimated, and a 'quasi-ideal' impulse response function with enhanced depth-resolution is achieved. We find out AR deconvolution is very suitable to deal with THz TDS signals, as the THz spectrum provides a large bandwidth with high SNR to establish a high-order AR model. It is also important to note that the deconvolution result is sensitive to the selection of parameters and SNR when a high-order AR model is utilized. This method is verified numerically and experimentally with a one-layered polymer coating with thickness of 22.5 µm, which cannot be resolved by FWDD. The thickness is successfully resolved by AR deconvolution, demonstrating that AR deconvolution enables us to beat the limitation of the THz wavelength and enhance the depth-resolution for resolving the optically thin layers in the THz regime.

The authors gratefully acknowledge the financial support of the Conseil Régional du Grand Est of the Fonds Européen de Développement Régional (FEDER), and of the Institut Carnot ARTS.

Fig. 1 .Fig. 2 .

 12 Fig. 1. The assumed impulse response function h 0 and the simulated received THz signal r 0 . The inset shows the actual THz reference signal f 0 .

Fig. 3 .

 3 Fig. 3. Comparison of the simulated deconvolution results obtained by AR deconvolution with p = 44, FWDD and the assumed impulse response function. The inset shows the AIC values with different model orders.

Fig. 4 .

 4 Fig. 4. (a) The deconvolution results based on FWDD and AR deconvolution. The inset shows the received THz signal reflected from the poly-coated steel sample; (b) The estimated frequency spectrum (up to 5 THz) based on autogressive spectral extrapolation.