Philippe Chartier

Mohammed Lemou

Florian Méhats

Gilles Vilmart

Highly-oscillatory problems with time-dependent vanishing frequency

Keywords: highly-oscillatory problems, time-dependent vanishing frequency, asymptotic expansion, uniform accuracy AMS subject classification (2010): 74Q10, 65L20

In the analysis of highly-oscillatory evolution problems, it is commonly assumed that a single frequency is present and that it is either constant or, at least, bounded from below by a strictly positive constant uniformly in time. Allowing for the possibility that the frequency actually depends on time and vanishes at some instants introduces additional difficulties from both the asymptotic analysis and numerical simulation points of view. This work is a first step towards the resolution of these difficulties. In particular, we show that it is still possible in this situation to infer the asymptotic behaviour of the solution at the price of more intricate computations and we derive a second order uniformly accurate numerical method.

Introduction 1.Context

Highly-oscillatory evolution equations of the form

U ε (t) := d dt U ε (t) = 1 ε AU ε (t) + f U ε (t) , U ε (0) = U 0 , 0 ≤ t ≤ T,
where T is a strictly positive fixed time, independent of ε, and where the operator A is supposed to be diagonalizable and to have all its eigenvalues in iZ (equivalently exp(2πA) = I), have received considerable attention in the literature, from both the point of view of asymptotic analysis [Per69, SV85, HLW06, CMSS12, CMSS15, CLM17] and the point of view of numerical methods [START_REF] Calvo | Numerical stroboscopic averaging for ODEs and DAEs[END_REF][START_REF] Ph | Multi-revolution composition methods for highly oscillatory differential equations[END_REF][START_REF] Castella | Stroboscopic averaging for the nonlinear Schrödinger equation[END_REF]. However, allowing the parameter ε to take values in a whole interval of the form ]0, 1], prevents the use of numerical methods constructed for specific regimes. As a matter of fact, standard methods 1 from the litterature [START_REF] Hairer | Solving ordinary differential equations. I[END_REF][START_REF] Hairer | Solving ordinary differential equations[END_REF] typically have error bounds expressed as powers of the stepsize h of the form2 error ≤ C h p ε q , p > 0, q > 0, where p is the order of the method and q is equal to p or p -1: while suitable for the regime ε close to 1, they require formidable computational power for small values of ε. At the other end of the spectrum, methods based on averaging and designed for small values of ε (see for instance [START_REF] Ph | Higher-order averaging, formal series and numerical integration I: B-series[END_REF]) typically admit error bounds of the form error ≤ C(h p + ε q ), p > 0, q > 0, where p is the order of the method and q is the order of averaging: they thus encompass an incompressible error for larger values of ε. In contrast, uniformly accurate methods [CLM13, [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Ph | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF] are robust schemes that are able to deliver numerical approximations with an error (and at a cost) independent of the value of ε ∈]0, 1], error ≤ Ch p .

In this paper, our objective is to construct uniformly accurate methods for equations whose frequency of oscillation depends on time. More precisely, we consider systems of differential equations of the form

U ε (t) = γ(t) ε AU ε (t) + f U ε (t) ∈ R d , U ε (0) = U 0 ∈ R d , 0 ≤ t ≤ T, (1.1)
where A ∈ R d×d and where the function f is assumed to be sufficiently smooth. The parameter ε again lies in the whole interval (0, 1] and the real-valued function γ is assumed to be continuous on [0, +∞). Many semi-classical models for quantum dynamics also assume the form of highly oscillatory PDEs with a varying frequency (which, once discretized in space, obey equation (1.1)), e.g. quantum models for surface hopping [START_REF] Chai | A multiband semiclassical model for surface hopping quantum dynamics[END_REF], graphene models [START_REF] Morandi | Wigner model for quantum trans-port in graphene[END_REF], or quantum dynamics in periodic lattice [START_REF] Morandi | Multiband wigner-function formalism applied to the zener band transition in a semiconductor[END_REF]. In such applications, the frequency γ may depend on time (and sometimes also on U ε and measures the gap between different energy bands, while the parameter ε is nothing but the Planck constant. We emphasize that the case of a varying frequency with a positive lower bound has been studied in [START_REF] Crouseilles | Nonlinear geometric optics methodbased multi-scale numerical schemes for a class of highly oscillatory transport equations[END_REF] for surface hopping, in [START_REF] Crouseilles | A micro-macro method for a kinetic graphene model in one-space dimension[END_REF] for graphene, and also in [START_REF] Hairer | Long-term analysis of the Störmer-Verlet method for hamiltonian systems with a solution-dependent high frequency[END_REF] where the long-term preservation of adiabatic quantities is established in a situation where the right-hand-side of equation (1.1) is Hamiltonian. However, the case where the frequency may become small (e.g. of the order of ε) or even vanish is more delicate and requires special attention from both analysis and numerical points of view. This is the reason why the main novel assumption in this article is that the function γ vanishes at some instant t 0 , or more precisely, that there exists (a unique) t 0 ∈ [0, T ] such that γ(t 0 ) = 0.

Our goal is to investigate problem (1.1) under these new circumstances, from both the asymptotic analysis (when ε → 0) and the numerical approximation viewpoints. For the sake of simplicity in this introductory paper, we assume that γ(t) is of the form

3 ∃p ∈ N * , ∀t ≥ 0, γ(t) = (p + 1)(t -t 0 ) p .
We emphasize that this situation is not covered by the standard theory of averaging as considered e.g. in [Per69, SV85, HLW06, CMSS10, CMSS15, CLM17], and that recent numerical approaches [CCMSS11, CLM13, CCLM15, CLMV18] are ineffective. All techniques therein indeed rely fundamentally on the assumption that γ(t) ≥ γ 0 uniformly in time, for some constant γ 0 > 0, and cannot be transposed to the context under consideration here4 .

Formulation as a periodic non-autonomous problem and main results

Upon defining u ε (t) = exp -(t-t 0 ) p+1 ε A U ε (t), the original equation (1.1) may be rewritten

uε (t) = F (t -t 0 ) p+1 ε , u ε (t) , u ε (0) = u ε 0 := exp - (-t 0 ) p+1 ε A U 0 , (1.2)
where F (θ, u) = e -θA f (e θA u) is 2π-periodic w.r.t. θ and smooth in (θ, u). We make the following assumption:

Assumption 1.1. The function f is twice continuously differentiable on R d and there exists M > 0 such that for all 0 < ε ≤ 1, equation (1.2) with t 0 ∈ [0, T ] has a unique solution on [0, T ], bounded by M , uniformly w.r.t. ε.

In the sequel, C will denote a generic constant that only depends on t 0 and on the bounds of ∂ α 2 F , α = 0, 1, 2, 3, on the set {(θ, u), θ ∈ T, |u| ≤ 2M }, where T = [0, 2π]. The aim of this work is now twofold. On the one hand, we show that, under mild and standard assumptions, an averaged equation (for (1.2) of the form

∀t ∈ [0, T ], uε (t) = F (u ε (t)), u ε (0) = u ε 0 (1.3)
persists (in F , function F is averaged w.r.t. the time variable).5 More precisely, we have the following theorem (see the proof in Section 2.2), which can be refined with the next-order asymptotic term (see Section 2.3).

Theorem 1.2. Suppose that Assumption 1.1 is satisfied and consider the solutions u ε (t), u ε (t) of problems (1.2), (1.3), respectively, on the time interval [0, T ]. Then, there exists ε 0 > 0 such that for all ε ∈]0, ε 0 [, and all t ∈ [0, T ],

|u ε (t) -u ε (t)| ≤ Cε 1 p+1 . (1.4)
Note that the bound |u ε (t) -u ε (t)| ≤ Cε obtained in the classical case [START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF] of a constant frequency (that is to say in the case where γ(t) = 1 in equation (1.1)), is degraded to (1.4) for p ≥ 1. For p = 0, both estimates coincide.

On the other hand, we construct in the case p = 1 a second-order uniformly accurate scheme for the approximation of u ε , that is to say a method for which the error and the computational cost remain independent of the value of ε (for more details on uniformly accurate methods, refer for instance to [START_REF] Ph | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Ph | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF]).

ω = 1 2π 2π 0 ω(θ)dθ.

Averaging results

We introduce the following function Γ :

[0, T ] → [0, S] with S = (T -t 0 ) p+1 + t p+1 0 , Γ(t) := t 0 |γ(ξ)|dξ = t p+1 0 + µ t (t -t 0 ) p+1 , µ t = sign(t -t 0 ) p = ±1,
and notice right away that Γ is invertible with inverse Γ -1 : [0, S] → [0, T ] given by Γ -1 (s) = s

1 p+1 0 + sign(s -s 0 ) |s -s 0 | 1 p+1 , s 0 = t p+1 0 .
Figure 1: The functions Γ (in blue) and Γ -1 (in red) with t 0 = 1 and T = 2 for p = 1, 2, 5.

Let us now consider for s = Γ(t) the function v ε (s) = u ε (t), which, for s = s 0 , satisfies

d ds v ε (s) = 1 Γ ′ • Γ -1 (s) uε Γ -1 (s) = 1 (p + 1)|s -s 0 | p p+1 F µs s -s 0 ε , v ε (s) (2.1) with initial condition v ε (0) = v ε 0 := u ε 0 , µ s = 1 if (s -s 0 ) p ≥ 0 -1 otherwise and F ±1 (θ, u) := F (±θ, u).
As an immediate consequence of Assumption 1.1, equation (2.1) has a unique solution on [0, S], bounded by M uniformly in 0 < ε ≤ 1. In this section, our aim is to show that there exists an averaged model for (2.1) of the form ∀s ∈ [0, S], vε (s) = 1

(p + 1)|s -s 0 | p p+1 F (v ε (s)), v ε (0) = v ε 0 , (2.2)
and then construct the first term of the asymptotic expansion of v ε (see Section 2.3). Note that, despite the singularity at s = s 0 of the right-hand side of (2.2), its integral formulation clearly indicates the existence of a continuous solution on [0, S].

Preliminaries

Let us introduce the following 2π-periodic zero-average functions

G ±1 (θ, u) = θ 0 (F ±1 (σ, u) -F (u))dσ - s 0 (F ±1 (σ, u) -F (u))dσ , and
H ±1 (θ, u) = θ 0 G ±1 (σ, u)dσ - s 0 G ±1 (σ, u)dσ .
Note that

1 2π 2π 0 F 1 (σ, u)dσ = 1 2π 2π 0 F -1 (σ, u)dσ = 1 2π 2π 0 F (σ, u)dσ = F (u)
which is the reason why F appears in lieu of F ±1 in the definition of G ±1 . It is clear that these functions and their derivatives in u are uniformly bounded: for |u| ≤ 2M , v ∈ R d and s ∈ R, we have

|G ±1 (s, u)| + |H ±1 (s, u)| ≤ C, |∂ 2 G ±1 (s, u)v| + |∂ 2 H ±1 (s, u)v| ≤ C|v|, (2.3) |∂ 2 2 G ±1 (s, u)(v, v)| + |∂ 2 2 H ±1 (s, u)(v, v)| ≤ C|v| 2 , (2.4)
where we have denoted ∂ 2 the partial derivative with respect to the variable u. We eventually define the function

∀u ∈ R d , ∀s ∈ R + , Ω ±1 (s, u) = +∞ s 1 σ p p+1 (F ±1 (σ, u) -F (u))dσ.
(2.5)

The following two technical lemmas will be useful all along this article.

Lemma 2.1. The function Ω ±1 is well-defined for all s ∈ R + and u ∈ R d . Moreover, for all u satisfying |u| ≤ 2M , all s ≥ 0 and all v ∈ R d , we have the estimates

|Ω ±1 (s, u)| ≤ C, |∂ 2 Ω ±1 (s, u)v| ≤ C|v|, |∂ 2 2 Ω ±1 (s, u)(v, v)| ≤ C|v| 2 . (2.6)
Restricting to strictly positive values of s, i.e. s > 0, we have furthermore

|Ω ±1 (s, u)| ≤ C s p p+1 , |∂ 2 Ω ±1 (s, u)v| ≤ C|v| s p p+1

,

(2.7)

and

Ω ±1 (s, u) + G ±1 (s, u) s p p+1 ≤ C s 1+ p p+1 , ∂ 2 Ω ±1 (s, u)v + ∂ 2 G ±1 (s, u)v s p p+1 ≤ C|v| s 1+ p p+1 .
(2.8)

Proof. We only prove the results for Ω ±1 as their adaptation to ∂ 2 Ω ±1 and ∂ 2 2 Ω ±1 is immediate. An integration by parts yields

Ω ±1 (s, u) = - G ±1 (s, u) s p p+1 + p p + 1 +∞ s 1 σ 1+ p p+1 G ±1 (σ, u)dσ,
where, from (2.3), the last integral is convergent and bounded by C s p p+1

. This yields the well-posedness of Ω ±1 for all s > 0 and (2.7). We now simply remark that for all s ≥ 0

Ω ±1 (s, u) = 1 s 1 σ p p+1 (F ±1 (σ, u) -F (u))dσ + Ω ±1 (1, u).
This gives the well-posedness for s = 0 and (2.6) can be deduced from (2.7) written for s = 1. A second integration by parts then gives

Ω ±1 (s, u) = - G ±1 (s, u) s p p+1 - p p + 1 H ±1 (s, u) s 1+ p p+1 + p p + 1 1 + p p + 1 +∞ s 1 σ 2+ p p+1 H ±1 (σ, u)dσ.

Previous integral is bounded by

C s 1+ p p+1
owing to (2.3) and this yields (2.8).

Remark 2.2. Since 1+s s p p+1 ≤ 2 for s ≥ 1, estimates (2.6) and (2.7) also imply for instance that for all s ≥ 0,

|Ω ±1 (s, u)| ≤ C (1 + s) p p+1 and |∂ 2 Ω ±1 (s, u)v| ≤ C|v| (1 + s) p p+1
.

In order to state next result, we now define, for any function φ :

T × R d → R d and for 0 ≤ a ≤ b ≤ S, the integral I ε (a, b) = 1 p + 1 b a 1 |σ -s 0 | p p+1 φ |σ -s 0 | ε , v ε (σ) dσ (2.9)
where v ε is assumed to be the solution of equation (2.1).

Lemma 2.3. For a given p ∈ N * , consider two smooth functions φ, ψ :

T × R d → R d satisfying the estimates |ψ(σ, u)| ≤ C and φ(σ, u) + ψ(σ, u) (1 + σ) p p+1 ≤ C (1 + σ) 1+ p p+1 , (2.10
)

for all θ ∈ T and all |u| ≤ M . If p = 1, we have ∀b ∈ [0, s 0 ], I ε (0, b) = √ ε 2 log ε + s 0 -b s 0 + ε ψ (v ε (b)) + O( √ ε), (2.11) ∀b ∈ [s 0 , S], I ε (s 0 , b) = √ ε 2 log ε b -s 0 + ε ψ (v ε (s 0 )) + O( √ ε), (2.12)
where averages are taken w.r.t. the first variable. If p ≥ 2, we have the estimate

∀0 ≤ a ≤ b ≤ S, |I ε (a, b)| ≤ Cε 1 p+1 .
(2.13)

Proof. Consider 0 ≤ b ≤ s 0 . A change of variables allows to write I ε (0, b) as

I ε (0, b) = ε 1 p+1 p + 1 s 0 ε s 0 -b ε 1 σ p p+1 φ (σ, v ε (s 0 -εσ)) dσ. Now, we split (p + 1)ε -1 p+1 I ε (0, b) = J 2 + J 3 + J 4 -J 1 into the sum of the four terms J 2 = s 0 ε s 0 -b ε 1 (1 + σ) 2p p+1 - 1 (σ(1 + σ)) p p+1 ψ (v ε (s 0 -εσ)) dσ, J 3 = s 0 ε s 0 -b ε 1 (σ(1 + σ)) p p+1 ( ψ -ψ) (σ, v ε (s 0 -εσ)) dσ, J 4 = s 0 ε s 0 -b ε 1 σ p p+1 r(σ, v ε (s 0 -εσ))dσ, J 1 = s 0 ε s 0 -b ε 1 (1 + σ) 2p p+1 ψ (v ε (s 0 -εσ)) dσ,
where we have denoted r(σ, u) = φ(σ, u)

+ ψ(σ,u) (1+σ) p p+1
. Owing to assumption (2.10) and 1

(1 + σ) 2p p+1 - 1 σ p p+1 (1 + σ) p p+1 ∼ - p p + 1 1 σ 3p+1 p+1
, σ → +∞, integrals J 2 and J 4 are absolutely convergent and uniformly bounded w.r.t. ε. As for J 3 , we use the relation

- (ψ -ψ )(σ, v ε (s 0 -εσ)) σ p p+1 (1 + σ) p p+1 = d dσ (κ (σ, v ε (s 0 -εσ))) + ε 1 p+1 (p + 1)σ p p+1 (∂ 2 κ F -µ ) (σ, v ε (s 0 -εσ))
where we have taken equation (2.1) into account with µ s = µ = (-1) p and

κ(s, u) = +∞ s (ψ -ψ )(σ, u) σ p p+1 (1 + σ) p p+1
dσ, in order to write J 3 as

J 3 = κ s 0 ε , v ε (0) -κ s 0 -b ε , v ε (b) + ε 1 p+1 (p + 1) s 0 ε s 0 -b ε 1 σ p p+1 (∂ 2 κ F -µ ) (σ, v ε (s 0 -εσ)) dσ
from which we may prove that J 3 is bounded (note indeed that ∂ 2 κ F -µ is bounded). For p > 1 it is clear that J 1 is bounded owing to (2.10) and finally, that I ε (0, b) is bounded. The contribution of J 1 for p = 1 is more intricate and requires to be decomposed as follows

J 1 = s 0 ε s 0 -b ε 1 1 + σ ψ (v ε (b))dσ + s 0 ε s 0 -b ε 1 1 + σ ψ (v ε (s 0 -εσ)) -ψ (v ε (b)) dσ = log s 0 + ε ε + s 0 -b ψ (v ε (b)) + s 0 ε s 0 -b ε 1 1 + σ ( ψ (v ε (s 0 -εσ)) -ψ (v ε (b))) dσ.
To estimate the second term, we use (2.1) and

s 0 -εσ ≤ b ≤ s 0 to get [ ψ (v ε (τ ))] s 0 -εσ b ≤ b s 0 -εσ 1 2 √ s 0 -τ ( ∂ 2 ψ F µ ) ( τ -s 0 ε , v ε (τ ))dτ ≤ C √ εσ so that s 0 ε s 0 -b ε ( ψ (v ε (s 0 -εσ)) -ψ (v ε (b))) 1 + σ dσ ≤ C √ ε s 0 ε 0 √ σ (1 + σ) dσ ≤ C √ s 0 .
We finally obtain that

I ε (0, b) = √ ε 2 log ε + s 0 -b s 0 + ε ψ (v ε (b)) + O( √ ε).
Mutatis mutandis, a similar conclusion holds true for the case a = s 0 and b ≥ s 0 as can be seen by writing the new value of

J 1 as b-s 0 ε 0 ψ (v ε (s 0 )) + ψ (v ε (s 0 + εσ)) -ψ (v ε (s 0 )) 1 + σ = log 1 + b -s 0 ε ψ (v ε (s 0 )) + O(1).

The averaged model

We are now in position to state the first averaging estimate, from which Theorem 1.2 follows by considering the change of variable Γ.

Proposition 2.4. Let v ε be the solution of problem (2.1) on [0, S], under Assumption 1.1.

Then, for all 0 < ε < ε 0 where ε 0 depends only on bounds on the derivatives of F , the solution v ε of the averaged model (2.2) exists on [0, S] and one has

∀s ∈ [0, S], |v ε (s) -v ε (s)| ≤ C ε 1 p+1 . (2.14)
Proof. The integral formulation of equation (2.1) reads

v ε (s) = v ε 0 + 1 p + 1 s 0 1 |σ -s 0 | p p+1 F (v ε (σ))dσ + R ε (s), (2.15)
where (with

µ σ = sign(σ -s 0 ) p ) R ε (s) = 1 p + 1 s 0 1 |σ -s 0 | p p+1 F µσ σ -s 0 ε , v ε (σ) -F (v ε (σ)) dσ, (2.16)
which is well-defined for all s ∈ [0, S]. From (2.5) with ς = sign(σ -s 0 ), σ = s 0 , we have

d dσ Ω ν |σ -s 0 | ε , v ε (σ) = ς ε (∂ 1 Ω ν ) |σ -s 0 | ε , v ε (σ) + (∂ 2 Ω ν ) |σ -s 0 | ε , v ε (σ) vε (σ) = - ς ε 1 p+1 |σ -s 0 | p p+1 F ςν σ -s 0 ε , v ε (σ) -F (v ε (σ)) + |σ -s 0 | -p p+1 (p + 1) (∂ 2 Ω ν ) |σ -s 0 | ε , v ε (σ) F µσ σ -s 0 ε , v ε (σ) , that is to say, taking ν = ςµ σ 1 |σ -s 0 | p p+1 F µσ σ -s 0 ε , v ε (σ) -F (v ε (σ)) = -ςε 1 p+1 d dσ Ω ςµσ |σ -s 0 | ε , v ε (σ) (2.17) + ςε 1 p+1 (p + 1)|σ -s 0 | p p+1 ∂ 2 Ω ςµσ |σ -s 0 | ε , v ε (σ) F µσ σ -s 0 ε , v ε (σ) ,
where we have used (2.1). For σ ≤ s ≤ s 0 we have µ σ = (-1) p = µ s , ς = -1 and therefore

R ε (s) = ε 1 p+1 p + 1 Ω -µs s 0 -s ε , v ε (s) -Ω -µs s 0 ε , v ε 0 (2.18) - ε 1 p+1 (p + 1) 2 s 0 1 (s 0 -σ) p p+1 ∂ 2 Ω -µs s 0 -σ ε , v ε (σ) F -µs s 0 -σ ε , v ε (σ) dσ
a relation from which we may deduce, using (2.6) and Assumption 1.1, that |R ε (s)| ≤ Cε 1/(p+1) . In particular, |R ε (s 0 )| ≤ Cε 1/(p+1) . As for s ≥ s 0 , we have µ σ = ς = 1 and thus

R ε (s) =R ε (s 0 ) + ε 1 p+1 p + 1 Ω 1 (0, v ε (s 0 )) -Ω 1 s -s 0 ε , v ε (s) (2.19) + ε 1 p+1 (p + 1) 2 s s 0 1 (σ -s 0 ) p p+1 ∂ 2 Ω 1 σ -s 0 ε , v ε (σ) F 1 σ -s 0 ε , v ε (σ) dσ
and we may again conclude from (2.6) and Assumption 1.1 that |R ε (s)| ≤ Cε 1 p+1 for s 0 ≤ s ≤ S and eventually for all 0 ≤ s ≤ S. Finally, we have on the one hand,

v ε (s) = v ε 0 + 1 p + 1 s 0 1 |σ -s 0 | p p+1 F (v ε (σ))dσ + O(ε 1 p+1 ),
and on the other hand,

v ε (s) = v ε 0 + 1 p + 1 s 0 1 |σ -s 0 | p p+1 F (v ε (σ))dσ,
as long as the solution of (2.2) exists. Assumption 1.1 and a standard bootstrap argument based on the Gronwall lemma then enable to conclude.

Next term of the asymptotic expansion

This section now presents how the estimate of Proposition 2.4 (analogously Theorem 1.2) can be refined by introducing an additional term of higher order in ε, namely ε 2 p+1 , in the asymptotic expansion.

Proposition 2.5. Let µ = (-1) p , and δ p = 1 if p = 1, δ p = 0 otherwise. Under Assumption 1.1, if we consider the solutions vε and wε of the equation (2.2) respectively on [0, s 0 ] and [s 0 , S] and with the respective initial conditions

vε (0) = v ε 0 - ε 1 p + 1 Ω -µ s 0 ε , v ε 0 , (2.20)
wε (s 0 = vε (s 0 ) + ε

1 p+1 p + 1 Ω 1 (0, vε (s 0 )) + Ω -µ (0, vε (s 0 )) - δ p ε 4 log ε ε + s 0 ∂ 2 G F (v ε (s 0 )),
and ṽε the continuous function defined by the following expressions:

s ≤ s 0 , v ε (s) = vε (s) + ε 1 p+1 p + 1 Ω -µ s 0 -s ε , vε (s) - δ p ε 4 log ε + s 0 -s ε + s 0 ∂ 2 GF (v ε (s)), s 0 ≤ s, v ε (s) = wε (s) - ε 1 p+1 p + 1 Ω 1 s -s 0 ε , wε (s) + δ p ε 4 log ε + s -s 0 ε ∂ 2 GF ( wε (s 0 )) + β ε
where

β ε = ε 1 p+1 p + 1 Ω 1 (0, wε (s 0 )) - ε 1 p+1 p + 1 Ω 1 (0, vε (s 0 )) , then we have ∀s ∈ [0, S], |v ε (s) -v ε (s)| ≤ C ε 2 p+1 .
(2.21)

Remark 2.6. In classical averaging theory (i.e. for γ(t) ≡ 1 or equivalently for p = 0), the solution v ε (s) of (2.1) is obtained as the composition of three maps (see for instance [START_REF] Perko | Higher order averaging and related methods for perturbed periodic and quasi-periodic systems[END_REF] or [START_REF] Sanders | Averaging methods in nonlinear dynamical systems[END_REF]): (i) a change of variable of the form v ε 0 + εϕ ε 0 (v ε 0 ) applied at initial time, (ii) the flow map at time s of a smooth differential equation whose vector field is of the form F + ε F ε and (iii) a change of variable of the form v ε 0 + εϕ ε s (v ε 0 ) applied time s. The ε 1 p+1 and log terms in (2.20) and in ṽε and wε are the counterpart of ϕ ε 0 (v ε 0 ) and ϕ ε s (v ε 0 ) in this more intricate situation.

Proof. In order to refine estimates (2.18) and (2.19) of R ε (s) obtained in the proof of Proposition 2.4, we rewrite them as s ≤ s 0 : R ε (s) = ε

1 p+1 p + 1 Ω -µ s 0 -s ε , v ε (s) -Ω -µ s 0 ε , v ε 0 -I -µ (0, s) , (2.22) s ≥ s 0 : R ε (s) = R ε (s 0 ) + ε 1 p+1 p + 1 Ω 1 (0, v ε (s 0 )) -Ω 1 s -s 0 ε , v ε (s) + I 1 (s 0 , s) , (2.23)
where the expression of I ε ν coincides with ), then, using (2.6)-(2.7), one has

I ε in Lemma 2.3 for φ(σ, u) = ∂ 2 Ω ν F ν (σ, u) and ψ(σ, u) = ∂ 2 G ν F ν (σ, u). If
∀ν = ±1, Ω ν s ε , x -Ω ν s ε , x ≤ Cε 1 p+1
and owing to (2.14), estimates vε (0

) -v ε (0) = O(ε 1 p+1 ) and wε (s 0 ) -vε (s 0 ) = O(ε 1 p+1
), and the Gronwall lemma, it stems that

∀0 ≤ s ≤ s 0 , v ε (s) -vε (s) = O(ε 1 p+1 ) and ∀s 0 ≤ s ≤ S, wε (s) -v ε (s) = O(ε 1 p+1 ).
As a consequence, v ε (s) can be replaced by vε (s) in (2.22) and by wε (s) in (2.23), up to O(ε 2 p+1 )-terms.

Case p > 1: Lemma 2.3 shows that the terms ε ), we thus have for s ≤ s 0

v ε (s) = v ε 0 + 1 p + 1 s 0 F (v ε (σ)) |σ -s 0 | p p+1 dσ + ε 1 p+1 p + 1 Ω -µ s 0 -σ ε , vε (σ) σ=s σ=0 + O(ε 2 p+1 ),
that is to say, by denoting

V ε (s) = v ε (s) -ε 1 p+1 p+1 Ω -µ s 0 -s ε , vε (s) , the equation V ε (s) -V ε (0) = 1 p + 1 s 0 F V ε (σ) + (v ε (σ) -V ε (σ)) (s 0 -σ) p p+1 dσ + O(ε 2 p+1 ) = 1 p + 1 s 0 F V ε (σ) dσ + ∂ 2 F (V ε (σ)) (v ε (σ) -V ε (σ)) (s 0 -σ) p p+1 dσ + O(ε 2 p+1 ), = 1 p + 1 s 0 1 (s 0 -σ) p p+1 F V ε (σ) dσ + O(ε 2 p+1 ),
where we have used Remark 2.2 to get the bound

s 0 1 (s 0 -σ) p p+1 ∂ 2 F (V ε (σ))Ω -µ s 0 -σ ε , vε (σ) dσ ≤ Cε 1 p+1 +∞ 0 1 (σ(1 + σ)) p p+1 dσ. From V ε (0) -vε (0) = O(ε 2 p+1
) and equation (2.2), we obtain by the Gronwall lemma

∀s ≤ s 0 , |ṽ ε (s) -v ε (s)| = |V ε (s) -vε (s)| ≤ Cε 2 p+1 .
For s ≥ s 0 , we write

v ε (s) = v ε (s 0 ) + 1 p + 1 s s 0 F (v ε (σ)) (σ -s 0 ) p p+1 dσ + (R ε (s) -R ε (s 0 )) = v ε (s 0 ) + 1 p + 1 s s 0 F (v ε (σ)) (σ -s 0 ) p p+1 dσ - ε 1 p+1 p + 1 Ω 1 σ -s 0 ε , wε (σ) σ=s σ=s 0 + O(ε 2 p+1 ). Denoting W ε (s) = v ε (s) + ε 1 p+1
p+1 Ω 1 s-s 0 ε , wε (s) , we have the simple equation

W ε (s) = W ε (s 0 ) + 1 p + 1 s s 0 F W ε (σ)) (σ -s 0 ) p p+1 dσ + O(ε 2 p+1 ),
and by comparing with equation (2.2), Gronwall lemma enables to conclude that W ε (s) -

wε (s) = O(ε 2 p+1 ) given that W ε (s 0 )-wε (s 0 ) = O(ε 2 p+1
) (by definition of wε (s 0 ) and W ε (s 0 ) and estimate (2.21) for s = s 0 ). The statement for s ≥ s 0 now follows from

β ε = O(ε 2 p+1 ).
Case p = 1: This case differs in that the terms

√ ε 2 I ε ν in (2.22
) and (2.23) are now of order ε log(ε) for s close to s 0 . This yields for s ≤ s 0

v ε (s) = v ε 0 + 1 2 s 0 F (v ε (σ)) √ s 0 -σ dσ + √ ε 2 Ω -µ s 0 -σ ε , vε (σ) σ=s σ=0 - √ ε 2 I -µ (0, s) + O(ε),
that is to say, by denoting

V ε (s) = v ε (s) - √ ε 2 Ω -µ s 0 -s ε , vε (s) + ε 4 log ε + s 0 -s ε + s 0 ∂ 2 G F (v ε (s)),
the equation

V ε (s) = V ε 0 + s 0 F V ε (σ) 2 √ s 0 -σ dσ + s 0 ∂ 2 F (V ε (σ)) 2 √ s 0 -σ (v ε (σ) -V ε (σ))dσ + O(ε) = V ε 0 + s 0 F V ε (σ) 2 √ s 0 -σ dσ + √ ε 4 s 0 ∂ 2 F (V ε (σ)) √ s 0 -σ Ω -µ s 0 -σ ε , vε (σ) dσ - ε 8 s 0 log ε+s 0 -σ ε+s 0 √ s 0 -σ ∂ 2 F (V ε (σ)) ∂ 2 G F (v ε (σ)) + O(ε) = V ε 0 + s 0 F V ε (σ) 2 √ s 0 -σ dσ + O(ε),
where we have used Lemma 2.3 again now with φ(σ, u) = ∂ 2 F (u) Ω -µ (σ, u) and ψ(σ, u) = ∂ 2 F (u) G -µ (σ, u), and noticed that ψ = ∂ 2 F G -µ = 0, to get rid of the second term of the second line. The third term may be bounded through an integartion by parts. We finally conclude by Gronwall lemma. For s ≥ s 0 , we get

v ε (s) = v ε (s 0 ) + 1 2 s s 0 F (v ε (σ)) √ σ -s 0 dσ - √ ε 2 Ω 1 σ -s 0 ε , wε (σ) σ=s σ=s 0 + √ ε 2 I 1 (s 0 , s) + O(ε),
that is to say, by denoting

W ε (s) = v ε (s) + √ ε 2 Ω 1 s -s 0 ε , wε (s) - ε 4 log ε ε + s -s 0 ∂ 2 G F ( wε (s 0 )),
the equation

W ε (s) = W ε (s 0 )+ s s 0 F W ε (σ) 2 √ σ -s 0 dσ- √ ε 4 s s 0 ∂ 2 F (W ε (σ)) √ σ -s 0 Ω 1 σ -s 0 ε , wε (σ) dσ + ε 8 s s 0 log ε ε+s-s 0 √ σ -s 0 ∂ 2 F (W ε (σ)) ∂ 2 G F ( wε (s 0 )) + O(ε) = W ε (s 0 ) + s s 0 F W ε (σ) 2 √ σ -s 0 dσ + O(ε),
where we have used equation (2.12) of Lemma 2.3, and we may conclude as before.

Corollary 2.7. Let µ = (-1) p , δ p = 1 if p = 1, δ p = 0 otherwise and τ 0 = t p+1 0 ε . Under Assumption 1.1, consider ūε 1 and ūε 2 , the solutions of

uε (t) = F (ū ε (t)), (2.24)
respectively on [0, t 0 ] and [t 0 , T ] with respective initial conditions

ūε 1 (0) = u 0 - ε 1 p+1 p + 1 Ω -µ (τ 0 , u ε 0 )
and

ūε 2 (t 0 ) = ūε 1 (t 0 ) + ε 1 p+1 p + 1 Ω 1 (0, ūε 1 (t 0 )) + Ω -µ (0, ūε 1 (t 0 )) + δ p ε 4 log (1 + τ 0 ) ∂ 2 G F (ū ε 1 (t 0 )).

Then we have

∀t ∈ [0, T ], |u ε (t) -u ε (t)| ≤ C ε 2 p+1 (2.25)
where ũε is the continuous function defined on [0, T ] by the following expressions:

0 ≤ t ≤ t 0 : u ε (t) = ūε 1 (t) + ε 1 p+1 p + 1 Ω -µ (τ, ūε 1 (t)) - δ p ε 4 log 1 + τ 1 + τ 0 ∂ 2 GF (ū ε 1 (t)), t 0 ≤ t ≤ T : u ε (t) = ūε 2 (t) - ε 1 p+1 p + 1 Ω 1 (τ, ūε 2 (t)) + δ p ε 4 log (1 + τ ) ∂ 2 GF (ū ε 2 (t 0 )) + β ε , with τ = |t-t 0 | p+1 ε and β ε = ε 1 p+1 p+1 Ω 1 (0, ūε 2 (t 0 )) -ε 1 p+1
p+1 Ω 1 (0, ūε 1 (t 0 )).

A micro-macro method

In this section, we suggest a micro-macro decomposition, analogous to the one introduced in [CLM17] and elaborated from the asymptotic analysis of Section 2. In a second step, we propose a uniformly accurate numerical method derived from this decomposition.

The decomposition method

Let u ε (t) be the solution of (1.2) and let u ε (t) be the approximation defined in Corollary 2.7, and consider the defect function

∆ ε (t) = u ε (t) -u ε (t), for t ∈ [0, T ]. (3.1) Proposition 3.1. Assume that f is of class C 2 and consider the solution u ε (t) of (1.2) on [0, T ]. For p ≥ 1, the function ∆ ε (t) defined by (3.1) satisfies ∀t ∈ [0, T ], |∆ ε (t)| ≤ Cε 2 p+1 , (3.2) ∀t ∈ [0, t 0 [∪]t 0 , T ], ∆ε (t) ≤ Cε 1 p+1 and if p = 1 ∆ε (t) ≤ C. (3.3) Proof.
By construction, ũε is continuous on [0, T ] and estimate (3.2) is nothing but (2.25). However, its derivatives are not continuous at t 0 . Hereafter, it is enough to consider t in [0, t 0 [ as the same arguments can be repeated for values in ]t 0 , T ]. From the expression of

u ε (t) = ūε 1 (t) + ε 1 p+1 p + 1 Ω -µ (τ, ūε 1 (t)) - δ p 4 log 1 + τ 1 + τ 0 ∂ 2 GF (ū ε 1 (t)), τ = (t -t0) p+1 ε , it stems by definition of Ω (see (2.5)) that ∆ε (t) = F -µ (τ, u ε ) -F -µ (τ, ūε 1 ) - ε 1 p+1 p + 1 ∂ 2 Ω -µ (τ, ūε 1 ) F (ū ε 1 ) - δ p √ ε 2 √ τ 1 + τ ∂ 2 G F (ū ε 1 ) + δ p ε 4 log 1 + τ 1 + τ 0 d dt ∂ 2 G F (ū ε 1 ) , (3.4)
where we have omitted t in u ε (t) and ūε

1 (t). Since |ū ε 1 (t) -u ε (t)| ≤ Cε 1 p+1 on [0, t 0 ] (and |ū ε 2 (t) -u ε (t)| ≤ Cε 1
p+1 on [t 0 , T ]), we have from Prop. 2.4 and Eq. (2.6), the following estimates

|F -µ (τ, u ε ) -F -µ (τ, ūε 1 )| ≤ Cε 1 p+1
and ε

1 p+1 p + 1 ∂ 2 Ω -µ (τ, ūε 1 ) F (ū ε 1 ) ≤ Cε 1 p+1 .
Besides, 2 √ τ ≤ 1+τ , |ε log ε| ≤ √ ε, and the first estimate of (3.3) is thus proven. Assuming now that p = 1 and using again equations (1.2) and (2.2), a second derivation leads to

∆ε (t) = - 2 √ τ √ ε ∂ 1 F (τ, u ε ) -∂ 1 F (τ, ūε 1 ) + ∂ 2 F (τ, u ε ) F (τ, u ε ) -2∂ 2 F (τ, ūε 1 ) F (ū ε 1 ) + ∂ 2 F (ū ε 1 ) F (ū ε 1 ) - √ ε 2 ∂ 2 2 Ω 1 (τ, ūε 1 ) F (ū ε 1 ), F (ū ε 1 ) - √ ε 2 ∂ 2 Ω 1 (τ, ūε 1 ) ∂ 2 F (ū ε 1 ) F (ū ε 1 ) + 1 -τ 2(1 + τ ) 2 ∂ 2 GF (ū ε 1 ) - √ ε 2 √ τ 1 + τ ∂ 2 G F (ū ε 1 ) + ε 4 log 1 + τ 1 + τ 0 d 2 dt 2 ∂ 2 G F (ū ε 1 ) .
Thanks to Assumption 1.1, Lemma 2.1 and (2.2), all the terms are clearly uniformly bounded, except the critical one in the first line, which requires more attention. We get

√ τ √ ε ∂ 1 F (τ, u ε ) -∂ 1 F (τ, ūε 1 ) ≤ C |t -t 0 | ε |u ε -ūε 1 | ≤ C |t -t 0 | ε |ū ε 1 -u ε | + C,
where we have used the result of Proposition 2.5, i.e. |u ε -u ε | ≤ Cε. It remains, using the expression of u ε , to observe that for t = t 0 , 0 < τ ≤ τ 0 so that owing to (2.7), we obtain

√ τ √ ε | u ε -ūε 1 | ≤ √ τ √ ε √ ε 2 |Ω 1 (τ, ūε 1 )| + ε 4 log 1 + τ 1 + τ 0 | ∂ 2 GF (ū ε 1 (t 0 ))| ≤ C √ τ √ ε √ ε √ τ + C τ τ 0 log 1 + τ 1 + τ 0 ≤ C.
This completes the proof.

A uniformly accurate first order numerical method

We are now in position to introduce uniformly accurate numerical schemes for (1.2). In this Section, we derive a uniformly accurate first-order method for p ≥ 1. Consider 0 = t 0 < . . . < t k < • • • < t N = T a subdivision of the interval [0, T ] containing the singularity t 0 , with h = max k=1,...,N (t k -t k-1 ). Inspired by the integral schemes in [START_REF] Ph | A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations[END_REF], we introduce the following method,

u k+1 = u k + t k+1 t k F (t -t 0 ) p+1 ε , u k dt (3.5)
Thanks to estimate (3.2) and the first estimate of (3.3), we obtain the following proposition.

3.2.

Assume that f is of class C 1 . Consider the solution u ε (t) of (1.2) on [0, T ], and the numerical scheme u k defined in (3.5). Then u k yields a uniformly accurate approximation of order one of the solution u ε (t k ). Precisely, there exist ε 0 > 0 and h 0 > 0 such that for all ε ≤ ε 0 and all h ≤ h 0 ,

|u k -u ε (t k )| ≤ Ch
for all t k ≤ T and where C is independent of ε and h.

The method (3.5) can be efficiently implemented numerically by using the Fourier expansion of the vector field F (θ, u), F (τ, u) = ℓ∈Z e iℓτ F ℓ (u).

The induction (3.5) then reads

u k+1 = u k + (t k+1 -t k )F 0 (u k ) + ℓ =0 ε ℓ 1/(p+1) F ℓ (u k ) t k+1 t k e iε -1 k(t-t 0 ) p+1 dt.
Using the change of variables s = ε -1 k(t -t 0 ) and introducing the notation Λ p (t) = +∞ t e is p+1 ds, we obtain the method (3.5) can be implemented numerically as

(3.6) u k+1 = u k + (t k+1 -t k )F 0 (u k ) + ℓ =0 ε ℓ 1/(p+1) F ℓ (u k ) Λ p ℓ ε 1/(p+1) (t k -t 0 ) -Λ p ℓ ε 1/(p+1) (t k+1 -t 0 ) .
Observe that the function Λ p (t) can be evaluated using the incomplete complex Gamma function Γ(ν, z) = +∞ z t ν-1 e -t dt where ν = 1/(p+1) for which efficient numerical packages exist.

A uniformly accurate second order numerical method

In this section, we introduce a scheme of uniform order two. The new method provides approximations (ū k , ∆ k ) of the pair (ū ε (t k ), ∆ ε (t k )). Assume that t 0 is one of the discretization points, i.e. t 0 = t k 0 for some k 0 . An approximation u k of u ε (t k ) is then derived by assembling the approximation u k of u ε (t k ) from formulas in Corollary 2.7 and eventually by setting u k = u k + ∆ k . Given that problem (2.24) is nonstiff, any second-order numerical scheme is suitable for the computation of ūk and thus of u k , and we simply choose here the Heun method ūk+1 = ūk + h 2

F (ū k ) + h 2 F ūk + h F (ū k ) .
As a consequence, we limit ourselves to the scheme for ∆ ε . Starting from

∆ ε (t k+1 ) = ∆ ε (t k ) + t k+1 t k F (τ (ξ), u ε (ξ) + ∆ ε (ξ)) dξ -( u ε (t k+1 ) -u ε (t k )), (3.7)
where τ (ξ) = |ξ-t 0 | p+1 ε , we consider at time t k+1/2 = t k +t k+1 2 the approximation

∆ k+ 1 2 = ∆ k + t k+ 1 2 t k F τ (ξ), u k + ∆ k dξ -( u k+ 1 2 -u k ).
Since the function u ε + ∆ ε = u ε has a bounded first time-derivative, the error associated to this scheme is of order O(h 2 ). Expanding F in Fourier series, we see that the scheme necessitates the computation of integrals of terms of the form e iℓξ 2 which may be easily computed numerically using the complex erf function. Now, for k < k 0 and t ≤ t 0 , we identify the smooth part of u ε (t) as

a ε (t) = ūε (t) + ∆ ε (t) - ε 4 log 1 + τ 1 + τ 0 ∂ 2 GF (ū ε (t)), so that u ε (t) = u ε (t) + ∆ ε (t) = a ε (t) + √ ε 2 Ω 1 (τ (t), ūε (t))
and, by Proposition 3.1 and its proof, it is clear that the second time-derivative of a ε is uniformly bounded. In order to approximate (3.7), we remark that

a ε (ξ) = a k + a k+1/2 -a k t k+1/2 -t k ξ -t k + O(h 2 ),
where setting ūk+1/2 = ūk + h 2 F (ū k ), we define for τ k+1/2 = τ (t k+1/2 ),

a k+1/2 = ūk+1/2 + ∆ k+1/2 - ε 4 log 1 + τ k+1/2 1 + τ 0 ∂ 2 GF (ū k+1/2 ).
Moreover, we have

∀(s, ŝ) ∈ R 2 + , Ω 1 (s, ūk ) -Ω 1 (ŝ, ūk ) = ŝ s F (σ, ūk ) -F (ū k ) √ σ dσ ≤ C| √ ŝ - √ s| so that √ ε and √ ε 2 Ω 1 (τ (ξ), ūε (ξ)) = √ ε 2 Ω 1 τ (ξ), ūk + √ ε 2 ξ -t k ∂ 2 Ω 1 τ (ξ), ūk F (ū k ) + O(h 2 ).
Therefore, denoting

b k = a k + √ ε 2 Ω 1 τ (t k ), ūk ,
our numerical scheme takes the form

∆ k+1 =∆ k + t k+1 t k F τ (ξ), b k dξ + t k+1 t k ξ -t k ∂ 2 F τ (ξ), b k a k+1/2 -a k t k+1/2 -t k dξ + t k+1 t k √ ε 2 ξ -t k ∂ 2 F τ (ξ), b k ∂ 2 Ω 1 τ (ξ), ūk F (ū k )dξ + t k+1 t k √ ε 2 ∂ 2 F τ (ξ), b k Ω 1 τ (ξ), ūk -Ω 1 τ (t k ), ūk dξ + u ε (t k+1 ) -u ε (t k ),
and has a truncation error of size O(h 3 ), uniformly in ε. As for k ≥ k 0 , we have

a k = ūk + ∆ k + ε 4 log 1 + τ 1 + τ 0 ∂ 2 GF (ū ε (t 0 )) + β ε , b k = a k - ε 1/2 2 Ω 1 τ (t k ), ūk , and
∆ k+1 =∆ k + t k+1 t k F τ (ξ), b k dξ + t k+1 t k ξ -t k ∂ 2 F τ (ξ), b k a k+1/2 -a k t k+1/2 -t k dξ - t k+1 t k √ ε 2 ξ -t k ∂ 2 F τ (ξ), b k ∂ 2 Ω 1 τ (t k ), ūk F (ū k )dξ - t k+1 t k √ ε 2 ∂ 2 F Γ(τ ) ε , b k Ω 1 τ (ξ), ūk -Ω 1 τ (t k ), ūk dτ + u ε (t k+1 ) -u ε (t k ).
According to the above computations, the uniform accuracy with second order of the proposed scheme may now be stated:

Proposition 3.3. Assume that f is of class C 2 Consider the solution u ε (t) of (1.2) on [0, T ], and the numerical scheme ( u k , ∆ k ) defined above. Then u k = u k + ∆ k yields a uniformly accurate approximation of order two of the solution u ε (t k ). Precisely, there exist ε 0 > 0 and h 0 > 0 such that for all ε ≤ ε 0 and all h ≤ h 0 ,

|u k -u ε (t k )| ≤ Ch 2
for all t k ≤ T and where C is independent of ε and h.

Numerical experiments

We test our method on the Hénon-Heiles system with solution U ε = (q 1 , q 2 , p 1 , p 2 ), U ε (t) = γ(t) ε p 1 , p 2 , -γ(t) ε q 1 -2q 1 q 2 , -q 2 -q 2 1 + q 2 2 , U ε (0) = (0.9, 0.6, 0.8, 0.5), with a time-varying parameter γ(t) = (p + 1)(t -t 0 ) p where t 0 is a zero of multiplicity p. The associated filtered system, satisfied by the variable u ε (t) ∈ R 4 defined by u ε (t) = (cos(θ)q 1 (t) -sin(θ) p 1 (t), q 2 (t), sin(θ)q 1 (t) + cos(θ)p 1 (t), p 2 (t)) , with θ = (t-t 0 ) 2 ε , takes the form (1.2) with F 1 (θ, u) = 2 sin θ u 1 cos θ + u 3 sin θ u 2 , F 2 (θ, u) = u 4 , F 3 (θ, u) = -2 cos θ u 1 cos θ + u 3 sin θ u 2 , F 4 (θ, u) = -u 1 cos θ + u 3 sin θ 2 +u 2 2 -u 2 .

We consider a time interval of length T = 1 and take t 0 = 1/3 as time where the oscillatory frequency vanishes. The reference solution is obtained using the matlab ode45 routine with a tiny tolerance. In Figures 2 and3, we have represented the error versus the stepsize of the numerical solution u k in (3.5) (uniform order 1) in cases where γ(t) has multiplicities p = 1 and p = 2 respectively. In Figure 4, we consider the method of Section 3.3 (uniform order 2) for p = 1. On the left pictures, the error is plotted as a function of the stepsize h, for fixed values ε ∈ {2 -k , k = 0, • • • , 11} where lines of slope 1 (Fig. 2 and 3) and slope 2 (Fig. 4) can be observed. On the right pictures, the error is plotted as a function of ε, for fixed values h ∈ {0.1/2 k , k = 0, • • • , 9}, which illustrates the uniform accuracy of the schemes with respect to ε. All curves are in perfect agreement with Propositions 3.2 and 3.3.

Conclusion

In this work, we have derived the first terms of the asymptotic expansion in ε of the exact solution of equation (1.2). As compared to standard averaging where γ is assumed to be bounded from below by a strictly positive constant, convergence towards the so-called averaged model is severally deteriorated for large values of p. For p = 1, the next term in the asymptotic expansion behaves quite unexpectedly as ε log(ε) when ε goes to zero and this seems to be the first time such a behaviour is revealed. Based on this asymptotic expansion, we have shown that it is possible to construct uniformly accurate numerical schemes of orders 1 for all p ≥ 1 and 2 for p = 1. Whether one may envisage to construct a uniformly accurate second-order method for p > 1 remains an open question and will be the subject of further investigations.

  x and x differ by an O(ε 1 p+1

Figure 2 :

 2 Figure 2: Method (3.5) (uniform order 1) for multiplicity p = 1. Error as a function of h for ε ∈ {2 -k , k = 0, • • • , 11} (left) and error as a function of ε for h ∈ {0.1/2 k , k = 0, • • • , 9} (right).

Figure 3 :

 3 Figure 3: Method (3.5) (uniform order 1) for multiplicity p = 2. Error as a function of h for ε ∈ {2 -k , k = 0, • • • , 11} (left) and error as a function of ε for h ∈ {0.1/2 k , k = 0, • • • , 9} (right).

Figure 4 :

 4 Figure 4: Method with uniform order 2 for multiplicity p = 1. Error as a function of h for ε ∈ {2 -k , k = 0, • • • , 11} (left) and error as a function of ε for h ∈ {0.1/2 k , k = 0, • • • , 9} (right).

Univ Rennes, INRIA, CNRS, IRMAR -UMR 6625, F-35000 Rennes, France. Philippe.Chartier@inria.fr

Univ Rennes, CNRS, IRMAR -UMR 6625, F-35000 Rennes, France. Mohammed.Lemou@univ-rennes1.fr

3 Univ Rennes, CNRS, IRMAR -UMR 6625, F-35000 Rennes, France. Florian.Mehats@univ-rennes1.

fr 4 Université de Genève, Section de mathématiques, 2-4 rue du Lièvre, CP 64, CH-1211 Genève 4, Switzerland, Gilles.Vilmart@unige.ch 1 Such as, for instance, the Runge-Kutta method used in the Matlab routine ODE45 (see the "Numerical experiments" Section 3.4).

The constant C here is independent of ε and h.

Note that applying an analytic time-transformation to (1.1) allows to consider more general analytic functions γ(t) and our analysis is not restricted to the polynomial case.

As a related recent work, we also mention the study[START_REF] Arnold | Stationary Schrödinger equation in the semiclassical limit: WKB-based scheme coupled to a turning point[END_REF] for the uniformly accurate approximation of the stationary Schrödinger equation in the presence of turning points which are spatial points used in quantum tunnelling models and where the spatial oscillatory frequency vanishes (analogously to our assumption γ(t0) = 0 with t0 = 0). However, the equation under consideration is linear and assumed to have an explicit solution on [t0, t1] for some t1 > 0 independent of ε. Beyond t1, the problem can be handled with a Wentzel-Kramers-Brillouin expansion, since the frequency is then lower bounded by positive constant.

Note that here as in the sequel, we denote the average of a function ω : T → R d by

Ω 1 τ (t k ), ūk ≤ Ch

Acknowledgements. The work of P.C., M.L., and F.M. is partially supported by the ANR project Moonrise ANR-14-CE23-0007-01, the INRIA associated team ANTIpODE and the ENS-Rennes project Muniq. The work of G.V. is partially supported by the Swiss National Science Foundation, grants No: 200020 178752 and 200021 162404.


Florentin Damiens

François Lott

Christophe Millet

Riwal Plougonven

Francois Lott

An Adiabatic Foehn Mechanism

Keywords: Foehn, mountain meteorology, mesoscale dynamics, gravity waves, critical level

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Foehn flows are characterized by downslope winds that are warm, dry and strong and that generally occur on the lee side of mountains [START_REF] Richner | Understanding and forecasting alpine foehn[END_REF]. For more than a century they have been studied worldwide (see for three different mountain ranges [Brinkmann(1974), [START_REF] Mcgowan | Observations of foehn onset in the southern alps, new zealand[END_REF], [START_REF] Drobinski | Foehn in the rhine valley during map: A review of its multiscale dynamics in complex valley geometry[END_REF]]), because they have strong economical impacts, for instance through increase in fire danger [START_REF] Takane | Formation mechanisms of the extreme high surface air temperature of 40.9c observed in the tokyo metropolitan area: Considerations of dynamic foehn and foehnlike wind[END_REF], [START_REF] Sharples | Foehn-like winds and elevated fire danger conditions is southeastern australia[END_REF], [START_REF] Cannon | Wrf simulation of downslope wind events in coastal santa barbara county[END_REF], ]), pollution [START_REF] Seibert | South foehn and ozone in the eastern alps -case study and climatological aspects[END_REF], ] or weather hazards. They can also actively contribute to the climate change through their local impacts on sea ice melting [START_REF] Elvidge | Foehn jets over the larsen c ice shelf, antarctica[END_REF], ].

One long-standing question about foehn concerns the origin of the warming of the air masses, and it is generally admitted that there are two dominant mechanisms that can explain it [START_REF] Richner | Understanding and forecasting alpine foehn[END_REF], ]. The first mechanism is along the line of meteorological textbooks [Seibert(2005), ] when the example of the foehn is used to illustrate the significance of moist processes by using thermodynamical diagrams. In this "diabatic" theory, the warming of the air masses results from the forced condensation that occurs upstream of the ridge, the large temperature and dry air result from the flow descent on the lee side. As observational evidence suggests that warming has a more adiabatic cause [Seibert(1990), ], a second mechanism, named "isentropic drawdown" is often proposed. In this mechanism, the foehn air comes from warmer and dryer altitudes upwind of the mountain because there is low level flow blocking. These two types of foehn, one diabatic with strong upstream ascent and one adiabatic with upstream flow blocking are sometimes referred to as Swiss and Austrian foehn respectively [START_REF] Wursch | Swiss and austrian foehn revisited: A lagragian based analysis[END_REF], ]. With the progress of Lagrangian analysis, the relative amplitude of these two mechanisms is now better established [START_REF] Smith | Orographic precipitation and air mass transformation: An alpine example[END_REF], [START_REF] Miltenberger | Revisiting the latent heating contribution to foehn warming: Lagrangian analysis of two foehn events over the swiss alps[END_REF], ] and it seems that the adiabatic mechanism is often more efficient that the diabatic one. The Lagrangian analyses also reveal that other irreversible processes can play a role, like the vertical mixing occurring when the air passes the mountain summit, and the radiative impact of the upstream clouds [START_REF] Elvidge | The causes of foehn warming in the lee of mountains[END_REF], ].

The fact that dynamics play such a central role in the foehn motivates the present paper. Actually, behind the apparently simple idea that upstream blocking elevates the air parcels before they abruptly descent on the lee side are hidden quite subtle nonlinear processes. These processes are controlled by the non dimensional mountain height, H N = HN U∞ , where H is the maximum mountain height, N the buoyancy frequency and U ∞ a scale for the incident wind. When H N is large, the vertical wavelength of the mountain waves U ∞ /N is small compared to the mountain height H, the mountain waves break near above the mountain top, and the resulting nonlinear dynamics produce strong downslope winds and foehn (see the review by [Durran(1990)] and recent simulations in [START_REF] Smith | Effects of inversion height and surface heat flux on downslope windstorms[END_REF]). The mixing produced by the breaking waves can induce a well mixed stagnant layer, which allows to establish a correspondence with the fully nonlinear hydraulic jump theory ([Baines (1998), Smith(1985)]). The presence of the stagnant layer also permits to argue that nonlinear reflections at the self-induced critical level present in the mixed layer can reflect downward the waves and amplify the response [START_REF] Laprise | On the structural characteristics of steady finite-amplitude mountain waves over bell-shaped topography[END_REF], [START_REF] Smith | Effects of inversion height and surface heat flux on downslope windstorms[END_REF]. Also, for large H N , part of the low level flow is blocked upstream making upstream blocking concomitant to the onset of foehn. All these efforts to interpret dynamically the onset of foehn in terms of breaking waves or blocking, hide the simple fact that some foehn can be produced by [Long(1953)]'s model for hydrostatic linear waves forced by a nonlinear boundary condition [START_REF] Lilly | The effects of terrain shape on nonlinear hydrostatic mountain waves[END_REF], Muraki(2011), ]. Interestingly, using such a model, [START_REF] Lilly | The effects of terrain shape on nonlinear hydrostatic mountain waves[END_REF]] predict that the wave response is enhanced for mountains with a small windward and large leeward slopes, which is a form of "isentropic drawdown" mechanism produced by blocking, if we assume that the blocked air can be represented dynamically by a gentle windward slope. Nevertheless, the results using [Long(1953)]'s type models are not much used to interpret foehn, essentially because when they produce foehn, they also predict intense wave breaking aloft which undermines all the linear interpretations of what occur downstream of the breaking zone (e.g which includes the foehn zone). As we shall see, low level wind shears can mitigate this upper level wave breaking because the mountain waves become less and less unstable with altitude when the wind increases.

Concerned with the recurrent problem in mountain meteorology theory that the near surface winds are ill defined ([Vergeiner(1971), Lott(2007)]) and can become very small, [Lott(2016)] (hereinafter L16) adapts the model of [Long(1953)] to incident shear flow that are null in z = 0. L16 shows that this induces a critical level dynamics near below the surface that produces downslope windstorm and foehn without upper level wave breaking. L16 also shows that downslope winds and foehn are favored when the surface flow is stable, i.e. when the surface Richardson number J = N 2 /U 2 z is large compared to its critical value, J = 0.25. As the model is linear, it also suggests that there is no need to interpret the onset of foehn in terms of upstream blocking, nonlinear internal wave reflections, or using hydraulic theory, which are all the ingredients that are often mentioned to explain adiabatic foehn.

The results in L16 have several limitations. The first is that the model in L16 is only adapted to situations when the shear layer depth, d, is much larger than the maximum mountain height, H, i.e. as long as the near critical level solutions in [START_REF] Booker | The critical layer for internal gravity waves in a shear flow[END_REF]] can be used to treat the critical level. The second is that L16 does not explain why the results found are so sensitive to J. The third is that it does not exclude the possibility that reflections can play a role on foehn, because linear reflections can be present. The fourth is that the inflow equations used in [Long(1953)] are only rigorously valid when the incident flow is uniform. The last is the representation of boundary layer absorptive effects: L16 uses Raighleigh drag and Newtownian cooling with a free slip boundary condition. Some earlier studies have explored more realistic boundary conditions, e.g. for turbulent flow over hills [Hunt et al.(1988a), Hunt et al.(1988b), ], or with viscous solutions [Lott(2007)]. In L16, these simplifications are made to ease the numerical treatment of the problem. They are in part justified by the fact that in the atmosphere the low level wind shears are often due to horizontal gradients in temperature, they do not solely result from boundary layer dynamics: they can extend well above it, and there the inviscid dynamics can be applied.

The purpose of the present paper is to fill these gaps by extending the L16's theory to taller mountains, and by validating it against fully nonlinear simulations. For this purpose, section 2 reformulates the results in L16 for a mountain of arbitrary height (not only for H << d) but still for small dissipation (when the linear drags correspond to a boundary layer depth z b << d). Section 3 exposes the dry foehn mechanism proposed and shows that it is entirely related to the near-surface critical level dynamics. Section 4 presents fully nonlinear simulations done with WRF and explains how the value of z b is chosen to enable a comparison of the theory with the numerical simulations. Section 5 characterizes systematically the results in terms of foehn intensity, downslope winds amplitude, and gravity wave stress. To check if the results remain robust beyond our simplified treatment of the boundary layer, it also presents tests where (i) the incident wind is not null at the surface and (ii) where a conventional atmospheric boundary layer scheme is used. Section 6 concludes and discusses further the relations with observations.

Theoretical model

To analyze the mountain waves produced by a stably stratified shear flow when the incident wind is null at z = 0, L16 considers the background flow profiles

U (z) = U ∞ tanh(z/d), N 2 (z) = const, (1)
incident on a 2-dimensional mountain, the height of which follows the Witch of Agnesi profile

h(x) = H 1 + x 2 2L 2 .
(2)

In (1), U (z) is the background horizontal wind, N (z) the Brunt Vaisala frequency, z the altitude, d is the vertical scale of the shear and U ∞ the incident wind maximum amplitude. In (2) H is the maximum mountain height, L its characteristic horizontal length, and x the horizontal coordinate. When scaling time by N -1 and distances by U ∞ /N the 2-dimensional non-rotating linear dynamics can be expressed in term of a non dimensional vertical velocity, w(x, z) of the form,

w(x, z) = +∞ -∞ f (k) ŵc (k, z)e ikx dk, (3)
where overbars denote dimensionless variables, and ŵc (k, z) is a canonical monochromatic solution of "unit" amplitude in the far field which satisfies the dissipative Taylor Goldstein equation,

d 2 ŵc dz 2 +    1 U -i z k √ J 2 - U zz U -i z k √ J -k 2    ŵc = 0. ( 4
)
In (3) the amplitude term f (k) is obtained by numerical inversion of the non-linear free-slip boundary condition:

w(x, h(x)) = U (h) + u(x, h(x)) ∂ x h(x), (5)
where u(x, z) is the horizonatl wind disturbance. In (4)-( 5)

U (z) = tanh z √ J , and h(x) = H N 1 + x 2 /2F 2 r (6) where J = N 2 d 2 U 2 ∞ , H N = HN U ∞ , and F r = LN U ∞ (7)
are the surface and minimum Richardson number, a non-dimensional mountain height, and a Froude number respectively. Still in (4) the dissipative vertical scale for each harmonics,

z k = z b 2 1 F r k + F r k . ( 8
)
has been expressed in terms of a global scale z b , with the first term in parenthesis resulting from Rayleigh drag and Newtonian cooling of coefficient z b /F r * √ J/2 and the second term from a viscous dissipation which acts in the x-direction only of coefficient z b F r / √ J/2 (see also the right hand sides of Eqs.4 in L16). They are introduced to regularize the critical level dynamics for all the harmonics, i.e. the longer and the shorter ones respectively. Note that z k is written differently than in L16 to make clear that when the formula is applied to the dominant wavenumber F r -1 , z k in (8) compares to the dissipative scale z b 1 . To construct ŵc , L16 uses for each k an exact solution of the inviscid version of (4), ŵinv (k, z), which is based on hypergeometric functions (see [START_REF] Lott | A transition from kelvin-helmholtz instabilities to propagating wave instabilities[END_REF]] and (34) in L16), and which asymptotic behaviors are

ŵinv (k, z >> 1) ≈ e -m z/ √ J , (9) ŵinv (k, z << 1) ≈ a 1 (k)z 1/2-iµ + a 2 (k)z 1/2+iµ (10)
where,

m = √ J |k 2 -1|, µ = |J - 1 4 |, (11)
and where a 1 (k) and a 2 (k) have analytical forms once the "unit" amplitude condition ( 9) is satisfied. Also, when k

2 < 1, m is changed in -isign(k)m
, where the sign is to ensure upward group speed. Near the surface, L16 also uses the asymptotic solution of the damped Taylor-Goldstein Equation ( 4) valid when z << 1:

ŵsrf (k, z) = a 1 (k)(z -z k ) 1/2-iµ + a 2 (k)(z -z k ) 1/2+iµ , (12)
where a 1 and a 2 are the same as in ( 10) to ensure matching. To establish how these aproximate solutions can be combined to provide a uniformly valid approximation of ŵc (k, z), Fig. 1 shows ŵinv (k, z) (thin solid), its approximation near z = 0 ((10), thick dashed) and ŵsrf (k, z) (thick dashed and dots). For the three values of J selected, J = 0.5, 2, 4, we see that below around z ≈ J/4 the inviscid exact solution ŵinv and its inviscid approximation on the right of (10) almost coincide whereas above this altitude the asymptotic damped approximation ( ŵsrf in ( 12)) and the asymptotic inviscid solutions (10) are very close. This suggests that we can take the altitude J/4 to make a transition between the exact inviscid solution and the asymptotic damped solution. In the following, this transition is made by patching these two functions using tanh tapers centered at z = J/4 and of depth ∆z = J/5, i.e. by writting,

ŵc (k, z) = 1 + tanh 5z J -5 4 2 ŵinv (k, z) + 1 -tanh 5z J -5 4 2 ŵsrf (k, z). ( 13
)
The thick gray curves in Figs. 1a,b,c show these new approximations, and illustrate how well they behave for all z, permitting a smooth passage from the damped solutions near the ground to the exact inviscid one in the far-field. Compared to L16, we use this approximation in the inversion of (5),

+∞ -∞ f (k) ŵc (k, h) -ûc (k, h) dh dx e ikx dk = U dh dx (14)
rather than using (12) to express ŵc and ûc in the terms between brackets in ( 14). This formally permits to consider more elevated ridges than in L16.

A dry foehn mechanism

As the numerical results obtained with this new formulation of the model are consistent with those in L16 we will only emphasize on the descents of warm air that characterize foehn. In the following we take a domain of length 800 spanned by 1024 equally spaced points, the Froude number F r = 5 and the boundary layer scale z b = 0.05. The three panels in the left column of Fig. 2 shows that when J increases the foehn effect increases, as expected from the results in L16, where the same behavior was discussed in terms of downslope windstorms. The reader is referred to this paper for more complete sensitivity tests to H N and z b , but it is worthwhile to note here that when J passes 1 typically, an isentropic surface initially at the altitude of the ridge is not much displaced upward when it approaches the ridge but easily reaches the surface along the downslope side of the ridge (lowest dashed line in Fig. 2c). This "foehn" effect becomes extremely pronounced when J = 4 in Fig. 2e: almost all the air located between H N < z < 2H N upstream of the ridge descent below the mountain top z < H N in the lee-side.

Near the surface downstream, the theory predicts highly distorted isentropes and potential temperatures that are below the minimum value upstream: this is not realistic so they are not shown but this indicates where the dynamics may trigger breaking.

If the foehn effect proposed is purely related to the critical level dynamics, then it should be present in a theory that is even simpler than the one used to produce the panels on the left in Fig. 2. For this purpose we next derive the theory when the wind shear U z and buoyancy frequency are both constant, in the hydrostatic approximation and when the boundary layer depth is identical for all harmonics z k = z b . As we will see in this "constant shear" case we can also formally eliminate the reflected waves, proving that they also play a minor role. To make clear that the dynamics in this case is only controlled by J and z b , the two parameters that control the critical level dynamics, we introduce a scaling suited to the constant wind shear configuration by taking for the horizontal and vertical scales those of the mountain. We also scale the horizontal velocity disturbance by the background velocity at the mountain top,

x = Lx, z = H z, u = U z H ũ, w = U z H 2 /L w, p = ρ r U 2 z H 2 p, and b = U 2 z H b. ( 15
)
In this new set-up and scaling, the Boussinesq dynamics writes

z∂ x ũ + w = -∂ x p, ∂ z p = b, z∂ xb + J w = 0, (16a)
∂ x ũ + ∂ z w = 0, w(z = h) = h + ũ ∂ h ∂ x , (16b) h
= 1 1 + x2 /2 , ( 16c
)
it is only controlled by the Richardson number J. In this setup, we know that for each harmonics the vertical velocity is exactly,

ŵ = z1/2+isign( k)µ , (17)
where we have only retain the upward propagating solution [START_REF] Booker | The critical layer for internal gravity waves in a shear flow[END_REF], ]: there is no reflected wave by construction. We then incorporate a dissipative scale zk = zb which we also take constant, a simplification that permits to define a canonical solution

ŵc (z, k) = z -isign( k)z b 1/2+isign( k)µ . ( 18
)
We will next use this solution to evaluate the wave field in (3) and the boundary condition in (5). To ease comparison, the results will be shown using the dimensionless variables of the more complete theory,

x = F r x, z = H N z, u = H N √ J ũ, w = H 2 N F r √ J w. ( 19
)
The results from the simplified theory are shown in the right panels of Fig. 2, the resolutions being identical to those used to build the corresponding panels on the left. To have comparable results nevertheless we have taken a boundary layer scale that is smaller than that used in Section 2, z b = 0.02, a difference that is only related to our choice to take a constant boundary layer scale (if we adopt the formula (8) we need to return to the value z b = 0.05 to obtain the same results). With this scale, the simplified theory gives almost exactly the same results as the more complete one, at least at low level. As we emphasized previously, in this simplified theory, there is no reflected waves by construction, proving that the onset of foehn can be entirely due to the near surface critical level dynamics i.e. a mechanism summarized mathematically by ( 18), a mechanism that exclude pre-conditionning of a well mixed layer by upper level wave breaking, non-linear and now linear wave reflections.

Finally to understand further why the near-surface critical level dynamics makes the foehn so sensitive to J it is worthwhile to notice that for large J, the [START_REF] Booker | The critical layer for internal gravity waves in a shear flow[END_REF]]'s solution ( 17) has an interpretation in terms of vertical wavenumber. To make this transparent we approximate it near the mountain top by writing z = H N + δz, which yields

z 1/2+iµ = z 1/2 e iµ ln z ≈ z 1/2 H iµ N e i µ H N δz . (20)
Near the mountain top, the disturbance has a characteristic wavenumber m ≈ µ/H N . Now, we know that mountain flow dynamics becomes nonlinear when the disturbance's vertical wavelength m -1 becomes smaller than the mountain height H N . The argument often given is that in this case the linear theory predicts that vertical velocities change sign before the air parcels reach the mountain top. If we apply this condition here, we have that nonlinear effect become large when H N /µ < H N , which approximately simplifies to the condition J > 1, which is the condition that favors foehn in the theory here and the downslope winds in L16.

Another point we make is that foehn is not associated with upper level wave breaking, i.e. that some form of wave breaking can occur below the mountain top without occurring aloft. To illustrate this qualitatively we next use the same line of arguments as before. We compare the wave vertical wavelength m B and the vertical displacement η B produced by the wave around a "breaking" altitude H B , yet to be found, and by writing z = H B + δz. This yields

m b ≈ µ/H B and η B ≈ H -1/2+iµ B H 3/2-iµ N . ( 21
)
respectively. Qualitatively, breaking occurs when |η B | ≥ m -1 B , which gives,

H B ≤ µ 2/3 H N . ( 22
)
Wave breaking is confined to the lowest layers, and inhibited above the mountain summit H N at least when µ < 1. This is the basic reason due to which the foehn mechanism described here is not associated with upper level wave breaking. Of course this effect simply translates the fact that in presence of shear, the vertical wavelength increases with altitude.

The result in ( 22) also illustrates that the depth of the layer over which breaking occurs increases with J. To a certain extent this suggests that the dissipative layer depth should increase with J. As we will see next this is in fact the case when we compare to fully nonlinear simulations.

Verification with fully nonlinear model

To extend the above results beyond the inflow-linear and Boussinesq approximations, we next use the WRF model in the 2D mountain flow configuration available on line2 [START_REF] Skamarock | A description of the advanced research wrf version 3[END_REF], ]. It is configured here for a 100km long and 12km height domain gridded by 1000 × 1000 points yielding resolutions of 100m and 12m in the horizontal and vertical directions respectively. At the lateral bounds the domain is open and below the model top a 6km deep absorbing layer prevents the gravity waves to be reflected downward. These quite high resolutions are to ensure that the near surface critical level dynamics is well represented. They require a time step ∆t = 0.5s. Note also that all the results from WRF shown next have been verified against simulations where the horizontal and vertical gridspacing are divided by two. As the WRF model is fully compressible, it is mandatory to recall that in it the stratification is represented by a Brunt Vaisala frequency given by N 2 = g θ 0 dθ 0 dz , where θ 0 is a background potential temperature. To ensure that it is constant and to minimize the vertical variations of the background fields that could result in partial reflections, we consider an isothermal atmosphere T r = 290K, yielding N 2 = κg 2 RTr = 3.10 -4 s -2 , where κ, R, and g have their conventional dry air and earth-values. We will then impose the background wind (1) and always take U ∞ = 10m/s, the minimum Richardson number will then be changed by varying the shear layer depth d. Also, the mountain height will be given by (2) with L = 2.8km to ensure F r = 5.

Finally, in all the WRF simulations presented the dissipative coefficients are set to 0 and the surface boundary condition is free slip. This does not mean that the simulations are inviscid, because irreversible and diabatic processes will necessarily occur if the model confirms the onset of convective overturnings like those found in the theoretical model in the previous section. All the experiments presented in the next section prove that the WRF dynamical core is stable and dissipative enough to handle these irreversible processes. When we return to the theory, it is a priori difficult to say how these irreversible processes can be translated in the theorety. Nevertheless, as the theory has a free parameter to represent dissipations, the boundary layer scale z b , we can vary it and try to identify a value for which there is a match between the nonlinear simulations and the theory. To sort out this issue we have undertaken a systematic comparison between the WRF model and the theoretical one, varying the boundary layer depth in the latter. We found a good agreement when we take a boundary layer depth varying like

z b / √ J = 0.05. ( 23
)
The 0.05 factor is purely empirical and little sensitive to the resolution in both models, but the √ J term suggests that between the simulations and the theory, the penetration of the boundary layer depth into the shear layer needs to be the same.

Fixed mountain height

To confirm that theory can match the fully nonlinear simulation, we present experiments with a fixed mountain height H = 250m, (H N = 0.5), and vary the minimum Richardson number between J = 0.5 and J = 5. Figure (3) shows the horizontal wind disturbances from the linear model (left) and the WRF model (right) for 3 different Richardson numbers (J = 0.5, 2, and 4). In each figure the contour interval is proportional to the background wind amplitude at the top of the Hill (≈ H N / √ J), a value representative of the disturbance amplitude needed to satisfy the lower boundary condition. The first striking result is that the left hand and right hand panels share strong similarities: the wave amplitude and phase line tilts with altitude are comparable. Also, when J increases, the phase line tilts with altitude become more and more vertical above the ridge, which means that the characteristic horizontal scale of the wave field decreases when J increases: this behavior is consistently reproduced by the theory and by the WRF model. In all panels, we also note that the horizontal wind disturbances are substantially larger on the downstream side (x > 0). This upstream/downstream difference becomes very pronounced near the ground, contrasting the two flanks of the mountain. It is this difference in horizontal winds that explains why low level shears can produce intense downslope winds (L16), and we see here that it is also present in the simulations.

Returning to the foehn effect, Fig. 4 shows the isentropic surfaces corresponding to the experiments shown in Fig. 3. Again, the results between the theory and WRF are quite consistent, some qualitative discrepancies emerge downstream where the simulations start to predict potential temperature jumps. These figures also illustrate another feature: the theoretical results shown in Fig. 2 and in Fig. 4 differ by the value of z b , which increases with J for the comparison with WRF (see eq. 23). Comparison of figures 2 and 4 shows that the increase of z b limits the onset of foehn.

Fixed penetration of the mountain into the shear

As said before, to match the WRF simulations, we had to keep constant in the theoretical model, the ratio between the shear layer depth √ J and the boundary layer scale z b . This limits the amplitude of foehn in the theoretical model, but does not make it disappear. To make the Foehn effect more evident in the simulations, we next keep the "penetration" of the mountain into the shear constant between the experiments by keeping H N / √ J constant. In this case, at the top of the hill, the incident wind is the same in all the experiments, and to a certain extent, it is only the stratification that is changed. Note also that in these experiments the ratio between the boundary layer depth and the mountain height z b /H N is also constant, which allows cleaner comparisons of results since z b now changes with √ J. The results are shown in Fig. 5 for (J = 0.5, H N = 0.25) and (J = 4, H N = 0.71), the corresponding plots for (J = 2, H N = 0.5) being already in Figs 4c-d. Again, a good match between the theory and WRF is found, which is important since we are now using nondimensional mountain height that are near 1. In both the theoretical model and in WRF, the foehn downstream is now slightly more substantial when J = 4 than when J = 0.5, but the difference are more in terms of horizontal distance over which the air descends on the lee-side. In the simulation with J = 4 the isentrope which is at the altitude of the mountain summit upstream rapidly descends below half the mountain top downstream, before returning to its upstream altitude after the passage through an hydraulic jump. Importantly both this behavior and the location of the jump are well predicted by the theory. When J = 0.5 the descent of air in WRF is less abrupt, in agreement with theory but the total descent is still about half the mountain height. Another point is worth noticing. In Figs. 5a-d, but this is also true in all the other Figures, there is very little uplift of the air parcels in front of the ridge, the most spectacular vertical displacements occur after the passage of the summit. This small uplift could suggest that diabatic effects should not much affect the results, at least when the atmosphere is not close to water vapor saturation.

In this paper we try to make the point that upper level wave breaking is not needed for foehn to occur. The basic argument is that the theory is linear, yet shows a very good agreement with the fully non linear simulations. Nonetheless, we can go one step further and diagnose if upper level wave breaking actually occurs in WRF and in the theory. The answer is simply that it does not,and to illustrate this the Figs. 5e-f display the isentropic surfaces over a much larger domain than in Figs. 5c-d. In no place, except in the downstream lee side, the isentropes become almost vertical. This is in contrast with what is found with the Long's model with uniform U (see Fig. 1f in L16). In it, downslope winds and foehn always occur in conjunction with upper level wave breaking.

5 More systematic comparisons and alternative configurations

Quantitative indices

To provide quantitive measures of our results we have built global indices characterizing our problem and conducted 64 pairs (theory versus WRF) of experiments, 32 corresponding to the constant height cases described in 4.1. and 32 to the constant penetration of the mountain into the shear cases described in 4.2. For the first index, we follow L16 and measure downslope windstorm intensity by,

A = M ax z<H N /2 0<x<5Fr u(x, z) U (H N / √ J) . (24)
The index A indicates to which extent the wind in the flow along the foothills (defined as the region with z < H N /2 and 0 < x < 5F r ) exceeds the background wind at the top of the hill (U (H/d)). In terms of foehn it is an important index because when the descent is rapid there is less time for diabatic effects to mitigate the warming of the air masses. The results in Fig. 6a somehow reproduce those in L16: the downslope wind intensity increases with stability J. However the increases are not as spectacular as in [Lott(2016)] where A much more easily reaches values above 1, but this is consistent with the fact that z b increases with J. This constraint is the price to pay for our theory to be consistent with the simulations. As our theory reproduces mountain waves quite well (see Fig 3), we have also evaluated the wave stress in z = 2 √ J (that is well above the maximum mountain height and the low level shear) and compare it to the linear hydrostatic wave drag produced by the mountain in (6) in a uniform flow of intensity U (H N ) = H N √ J :

F z = - ∞ -∞ u w dx π 4 H 3 N √ J . ( 25
)
Again, the curves in Fig. 6b show that the theory and the nonlinear simulations are in good agreement. Interestingly, when the penetration of the mountain into the shear stays constant and J increases, the normalized wave stress F z increases well beyond 1. This behavior is reminiscent of the high drag states often found in non linear mountain flow theory [START_REF] Bacmeister | On high drag states of nonlinear stratified flow over an obstacle[END_REF], ].

Finally, to characterize foehn, we identify the highest altitude from which are coming the isentropic surfaces that arrive along the downstream foothill, and normalize by the mountain maximum altitude. In our non-dimensional formalism where the background potential temperature varies as z exactly, this is conservatively given by,

F hn = M ax z<H N /2 0<x<5Fr θ + z H N . (26)
Interestingly, in our models this factor is always near 1. In other words the dynamics forces isentropic surfaces that are upstream at the altitude of the ridge to descend along the foothills downstream. Hereagain simulations and theory compare well. In both, the foehn intensity increases with J, but as for the downslope winds intensity, it is less pronounced than when z b is kept small and constant in the theory (not shown).

Alternative configurations

In the real atmosphere the low level shears have at least two causes: the thermal gradients due to the presence of fronts and the effect of the boundary layer. The first imposes shears that extend well above the boundary layer, the second forces the winds to vanish at the surface. There is no reason, except simplicity, to consider that the shears associated with these two processes have similar values. In this sense, we could say that our almost inviscid interpretation of the wave dynamics above the boundary layer is only adapted to the shears produced via thermal wind balance, Raighleigh friction and Newtownian cooling being an ad hoc parameterizations of the boundary layer absorptive properties, an absorption that becomes very efficient when the incident wind is null at z = 0. To show that our results have more generality than this apparently "ad-hoc" setup, we next propose three experiments. In the first two, we keep the same models and setup as in Figs. 5c and 5d but shift vertically the mountain height by H N /10, e.g. by writing,

h(x) = h = H N 1 + x 2 /2F 2 r + H N 10 . ( 27
)
Doing so there is no need to rewrite the theory but, from a dynamical point of view there is a major change: the incident wind is no longer exactly null at the surface. The results in Figs. 7a and7b show that in this case the theory and WRF simulations are still consistent and that the foehn effects are almost identical to those when when the incident wind is null at z = 0 in Figs. 5c and5d respectively.

In the third experiment, we return to the exact setup of WRF used in Fig. 5d, but instead of the free slip boundary layer condition we adopt a planetary boundary layer scheme that has been extensively used in WRF. The chosen configuration activates the vertical diffusion described in [START_REF] Hong | Nonlocal boundary layer vertical diffusion in a medium-range forecast model[END_REF]] with the surface conditions evaluated using similarity theory. At the surface we impose a constant potential temperature, with value equal to the upstream temperature value at z = 0 and take a roughness length z 0 = 10cm. The results for the potential temperature and after a non-dimensional time t ≈ 100, are shown in Fig. 7c. Again the results in term of foehn compare well with those in Fig. 5d, and this despite the fact that now the low level winds and temperature adapt themselves to more conventional surface boundary conditions.

Discussion

The central result of this paper is the resemblance between Fig. 5c and Fig. 5d. It shows that isentropic "drawdown", a dry mechanism of foehn, where air parcels do not ascend significantly on the upstream side of mountains before descending abruptly on the lee side can be produced in a model without inflow nonlinear effects. The mechanism at work is related to the near surface critical level dynamics that occurs for mountain waves, the critical level being located just below the lower boundary. In this case, the inflow linear dynamics anticipate well where the wave breaks as is often the case in gravity wave theory. Not surprisingly, in the presence of low level wind shear, mountain wave breaking is favored in regions located below the mountain top rather than above.

In essence, the results obtained illustrate the significance of low level wind shears and stability on mountain flow dynamics, something often noticed in mountain meteorology, at least in the context of trapped lee waves [START_REF] Reinecke | Initial-condition sensitivities and the predictability of downslope winds[END_REF], [START_REF] Georgelin | On the transfer of momentum by trapped lee waves. case of the iop3 of pyrex[END_REF]. In this context it is interesting to note that there are many field experiments where upstream soundings show low level shears [START_REF] Doyle | An intercomparison of t-rex mountain-wave simulations and implications for mesoscale predictability[END_REF]), Sheridan et al.(2007), ], generally related to advancing fronts [START_REF] Lothon | Experimental study of five fhn events during the mesoscale alpine programme: From synoptic scale to turbulence[END_REF], ], when the low level shear is associated with thermal wind balance. It is also worth noticing that there are cases where foehn is observed in the absence of upper level wave breaking [START_REF] Flamant | Gap flow in an alpine valley during a shallow south fhn event: Observations, numerical simulations and hydraulic analogue[END_REF]]. Other aspects of the theory also have reasonable behavior. According to [START_REF] Richner | Understanding and forecasting alpine foehn[END_REF]] a forecaster's rule says that foehn winds rarely descend more than 2000m below the mountain crest. In all figures when the foehn is strong (for instance when J ≥ 2 in Figs. 45) we note a comparable behavior, the foehn never extend down to z = 0 in the lee side in these cases. Also, downstream where the isentropic surface become vertical, there is always near the surface a pool of relatively cold air. It is interesting that such a pool has often been viewed as an active part of the foehn dynamics [START_REF] Lothon | Experimental study of five fhn events during the mesoscale alpine programme: From synoptic scale to turbulence[END_REF], ], we see here that it can be a consequence of it.

All these remarks are not to say that we want here to replace the existing theories of adiabatic foehn by a completely new one, the message is rather to complement them. Actually, our upstream low level shear can be viewed as a form of blocking, since it imposes slow incident low level flows. Also, upper level mountain wave breaking tend to decelerate the low level flow, again our background profiles with slow low level winds can be viewed as a result of this process. Finally, our results show that linear models forced by a nonlinear boundary condition potentially have some skill, which can turn out to be useful since simplified linear models are still sometimes used to predict mountain wave operationally at a cheap numerical cost [Vosper(2003), [START_REF] Sheridan | Mountain waves in high resolution forecast models: Automated diagnostics of wave severity and impact on surface winds[END_REF], ]. Finally, one may note that the equations solved in the numerical simulations and in the linear theory differ in their treatment of boundary layer dissipation. For consistency one may wish that the WRF model and the theory be adapted to use the same dissipative processes. Nevertheless, we explored in Section 5 the sensitivity of the results to the precise boundary condition at the surface (non-zero surface wind or including a boundary layer parameterization) and these tests suggested that the foehn effects are robust to

Figure 1 :

 1 Figure

Figure 2 :

 2 Figure 2: Isentropic surfaces predicted by the linear model when H N = 0.5, F r = 5 and for different value of J. Tanh incident wind and z b = 0.05: a) J=0.5, c) J=2, e) J=4. Constant shear and hydrostatic set-up in (16), z b = 0.02: b) J=0.5, d) J=2, e) J=4. In all panels the thick curve represent the incident wind U and the thick dashed curves the isentropes that initially are at the altitudes z = H N and z = 2H N .

Figure 3 :

 3 Figure 3: Non-dimensional Horizontal wind disturbances (u) predicted by the theory (a, c, e) and WRF (b, d, f), for J = 0.5, 2, and 4, respectively. In all simulations H N = 0.5 and in the linear model in (a), (c), and (e) the boundary layer scale varies as z b = 0.05 √ J. In each panel the contour interval is proportional to the background wind at the top of the hill, contours= 0.2 * H N / √ J.

Figure 4 :

 4 Figure4: Isentropic surfaces predicted by the theory and WRF, same parameters as in Fig.3

Figure 5 :

 5 Figure 5: Same as Fig. 4 but leaving H N / √ J =cte: a) and b) J = 0.5, H N = 0.25; c) and d) J = 4, H N = 0.71. e) and f) are enlarged views of c) and d).

Figure 6 :

 6 Figure 6: Global indices measuring a) the downslope wind amplitudes (24), b) the normalised wave stress (25), and c) the foehn intensity (26). The results are from 64 numerical experiments, 32 uses with WRF and 32 the "inflow linear" theory (LT). Half the simulations are done with constant height H N = 0.71, the other half with constant penetration of the mountain into the shear H N / √ J = 0.4.

Figure 7 :

 7 Figure 7: Sensitivity test: a) and b) are as Fig. 5c and 5d with the surface h shifted vertically by H N /10. c) is as Fig. 5d but with a fully parameterized boundary layer.

For consistency, note also that the boundary layer depth z B used in L16 is related to the one here by z B = 5z b .

www2.mmm.ucar.edu/wrf/users/downloads.html

Acknowledgements

This work was supported by the European Commission's project ARISE2 (Grant Agreement 653980) and the Comissariat à l'Energie Atomique (CEA).

such changes.


Jabrane Karkouri

Fabien Millioz

Thomas Troalen

Rémy Prost

Magalie Viallon

Hélène Ratiney

ACCELERATED SPIRAL CHEMICAL SHIFT IMAGING FOR PROTON DENSITY AND T2*FAT-WATER QUANTIFICATION

whether they are published or not. The documents may come

INTRODUCTION

Usual chemical shift encoding water-fat MRI methods require strong constraints and prior knowledge of the water-fat spectral model to obtain fat quantification and approximation of a global T2*. In contrast, traditional Magnetic Resonance Spectroscopic Imaging (MRSI) method uses less constrained spectral model but faces long acquisition time [START_REF] Adalsteinsson | Volumetric spectroscopic imaging with spiral-based k-space trajectories[END_REF][START_REF] Weis | Lipid content in the musculature of the lower leg: evaluation with high-resolution spectroscopic imaging[END_REF]. We propose an accelerated spiral MRSI method based on temporal irregular under-sampling [START_REF] Karkouri | Time samples selection in spiral acquisition for sparse magnetic resonance spectroscopic imaging[END_REF] and demonstrate its applicability in vivo on the calf.

2. THEORY AND METHOD

Design of the temporal undersampling

Obtaining a spectrum from irregularly undersampled time domain signal is possible with least square (LS) reconstruction and assuming a known spectral support but results in a noise amplification. We minimize this noise amplification using the sequential backward selection (SBS) algorithm for time sample selection.

SBS based spiral MRSI

In spiral MRSI, multiple k-spaces are acquired over time with spiral trajectories [START_REF] Karkouri | Time samples selection in spiral acquisition for sparse magnetic resonance spectroscopic imaging[END_REF]. Spiral trajectory determines the dwell time which can be too long compared to desired temporal resolution. As a result temporal interleaves are needed which lengthens the acquisition time. Even with this constraint SBS subsampling result in a shorter acquisition time. Spatially distributed spectra are obtained from the k-t data with a gridding and a 2D Fast Fourier Transform (FFT) for the spatial reconstruction, followed by a 1D FFT for the fully sampled case or LS minimization for the irregularly undersampled case. 3. DISCUSSION AND CONCLUSION Experiments showed proper spectra reconstructions in the case of our under-sampled method, as T2*, fat and water maps. This method would have a larger impact on time gain reduction both for multi-slice, larger matrices applications. Red arrow in Fig. 1 highlights physiologic elevated T2* in previously exercised muscle.

2 . 3 .

 23 Spectrum quantificationWater and fat amplitude (Aw, Af) and damping factor (dw and df) of were fitted on each time domain signal using a exponential decay model, enabling to create proton density (proportional to A) and T2* equal to (1/d) maps. 2.4. in vivo application: human calf A Siemens PRISMA 3T machine has been used, with a 1H surface coil, a matrix 32x32, 6 spatial interleaves, FOV 20cm, 2.5ms spiral duration, 255ms readout duration, TR 2s for an acquisition time of 1min12s for the fully temporal sampled spiral MRSI, 48s for the under-sampled (reduction factor 2). T2* values are ranged between [15, 45] ms for water, [15, 80] ms for fat and the mean fit error is 10 (compared to reconstructed signal).

Fig. 1 .Fig. 2 .

 12 Fig. 1. Fat, water (amplitude in u.a) and T2* (in ms) maps, in fully sampled (first row) and under-sampled (second row) time domain acquisitions