Youcef Amirat
email: youcef.amirat@uca.fr.

Arnaud Münch
email: arnaud.munch@uca.fr.

Internal layer intersecting the boundary of a domain in a singular advection-diffusion equation

Keywords: Asymptotic analysis, Singular perturbation, Internal and boundary layers, Sobolev estimates AMS subject classification (2020): 35C20, 35K67

We perform an asymptotic analysis with respect to the parameter ε > 0 of the solution of the scalar advection-diffusion equation y ε t +M (x, t)y ε x -εy ε xx = 0, (x, t) ∈ (0, 1)×(0, T ), supplemented with Dirichlet boundary conditions. For small values of ε, the solution y ε exhibits a boundary layer of size O(ε) in the neighborhood of x = 1 (assuming M > 0) and an internal layer of size O(ε 1/2 ) in the neighborhood of the characteristic starting from the point (0, 0). Assuming that these layers interact each other after a finite time T > 0 and using the method of matched asymptotic expansions, we construct an explicit approximation P ε satisfying y ε -P ε L ∞ (0,T ;L 2 (0,1)) = O(ε 1/2 ). We emphasize the additional difficulties with respect to the case M constant considered recently by the authors.

Introduction. Problem statement

Let T > 0 and Q T := (0, 1) × (0, T ). This work is concerned with the scalar advection-diffusion equation      y ε t (x, t) + M (x, t) y ε x (x, t) -εy ε xx (x, t) = 0, (x, t) ∈ Q T , y ε (0, t) = v(t), y ε (1, t) = 0, t ∈ (0, T ), y ε (x, 0) = y 0 (x), x ∈ (0, 1), [START_REF] Amirat | Asymptotic analysis of an advection-diffusion equation and application to boundary controllability[END_REF] where ε ∈ (0, 1) is the diffusion coefficient and M (x, t) > 0 is the transport velocity. For any initial data y 0 ∈ H -1 (0, 1) and Dirichlet condition v ∈ L 2 (0, T ), there exists a unique solution y ε ∈ L 2 (Q T ) ∩ C([0, T ]; H -1 (0, 1)). This apparently simple partial differential equation appears in many situation as it is a prototype of models where the diffusion coefficient is small compared to the others. As explained in [START_REF] Cheng | New approximation algorithms for a class of partial differential equations displaying boundary layer behavior[END_REF], this model can notably be seen as an embedded system of the Navier-Stokes system with non-characteristic boundary condition and viscosity coefficient equals to ε. It also can be seen as a regularization of a transport equation and for this reasons well employed in numerical analysis (see notably [START_REF] Deuring | L 2 -stability independent of diffusion for a finite elementfinite volume discretization of a linear convection-diffusion equation[END_REF][START_REF] Stynes | Uniformly convergent difference schemes for singularly perturbed parabolic diffusion-convection problems without turning points[END_REF] and the references therein) when one wants to obtain robust numerical approximation uniformly with respect to ε small. Equation [START_REF] Amirat | Asymptotic analysis of an advection-diffusion equation and application to boundary controllability[END_REF] appears in models of miscible displacement of compressible fluids in porous media, with small molecular diffusion and dispersion coefficients, see [START_REF] Bear | Introduction to modeling of transport phenomena in porous media[END_REF]. Last, it also appears in the context of exact boundary controllability when one wants to steer to zero the solution y ε with a uniform control v: we mention the seminal works [START_REF] Coron | Singular optimal control: a linear 1-D parabolic-hyperbolic example[END_REF][START_REF] Guerrero | Singular optimal control for a transport-diffusion equation[END_REF] and the recent paper [START_REF] Laurent | On uniform observability of gradient flows in the vanishing viscosity limit[END_REF].

We are interested in this work with a precise asymptotic description of the solution y ε when ε is small. This problem has been the subject of several studies in the last decades in the case for which the transport velocity is constant and the equation is defined over R + × (0, T ). We refer to [START_REF] Shih | A novel uniform expansion for a singularly perturbed parabolic problem with corner singularity[END_REF][START_REF] Shao | Asymptotic analysis and domain decomposition for a singularly perturbed reaction-convectiondiffusion system with shock-interior layer interactions[END_REF]. The case of a transport velocity depending only on the time variable has been formally discussed in [START_REF] Plaschko | Matched asymptotic approximations to solutions of a class of singular parabolic differential equations[END_REF] and deeper analyzed in [START_REF] Shih | On a class of singularly perturbed parabolic equations[END_REF].

For bounded domains with respect to the space variable, the asymptotic analysis is quite involved as several singular layers may appear in Q T and interact each other. The constant transport velocity case has been analyzed in [START_REF]Asymptotic analysis of an advection-diffusion equation involving interacting boundary and internal layers[END_REF]. More precisely, if we denote the characteristic t → X(t; x, s) through (x, s) ∈ Q T as the solution in Q T of dX dt = M (X, t), X(s; x, s) = x,

then the violation of the compatibility conditions between y 0 and v at the point (0, 0) for which y 0 (0) = v(0) induces a thin inner region (called internal layer) of size O(ε 1/2 ) in the vicinity of the characteristic {(x, t) ∈ Q T , x -X(t; 0, 0) = 0} where the solution y ε exhibits rapid variations. Thus, if this characteristic gets arbitrarily close to the line x = 1 in a finite time T 1 > 0 unique solution of the equation

X(T 1 ; 0, 0) = 1, (3)
then the internal layer interacts with the usual boundary layer of size O(ε) living along x = 1 and induced by the Dirichlet condition. Figure 1 provides a geometric description of this phenomenon in two cases: the case in Figure 1-Left for which the function M is constant, i.e. M (x, t) = M > 0 leading to a linear characteristic of equation x -M t = 0. The internal layer lives in the red zone {(x, t) ∈ Q T ; |x -M t| ≤ ε 1/2 } and intersects the boundary layer occupying the blue zone {(x, t) ∈ Q T , |x -1| ≤ ε} in a small neighborhood of the point (1, 1/M ) such that T 1 = 1/M . This case has been extensively studied in [START_REF]Asymptotic analysis of an advection-diffusion equation involving interacting boundary and internal layers[END_REF] (see also [START_REF] Amirat | Asymptotic analysis of an advection-diffusion equation and application to boundary controllability[END_REF][START_REF]On the controllability of an advection-diffusion equation with respect to the diffusion parameter: asymptotic analysis and numerical simulations[END_REF]). Figure 1-Right corresponds to more general situation for which the function M is not constant. The characteristic starting from the point (0, 0) has the equation x -X(t, 0, 0) = 0. In this paper, we extend the analysis given in [START_REF]Asymptotic analysis of an advection-diffusion equation involving interacting boundary and internal layers[END_REF] devoted to a constant velocity M (x, t) = M > 0 and perform an asymptotic analysis of the solution y ε of (1). Precisely, our main result reads as follows.

ε √ ε x t 0 0 1 T ≥ 1 M t = 1 M y ε (x, 0) = y 0 (x) y ε (0, t) = v(t) y ε (1, t) = 0 ε √ ε x -a ( t ) = 0 x t 0 0 1 T ≥ T 1 t = T 1 y ε (x, 0) = y 0 (x) y ε (0, t) = v(t) y ε (1, t) = 0
Theorem 1.1. Assume y 0 ∈ C 2 ([0, 1]), v ∈ C 2 ([0, T ]) and M ∈ C 2 (Q T , R + ). For ε > 0 small enough, let y ε be the solution of [START_REF] Amirat | Asymptotic analysis of an advection-diffusion equation and application to boundary controllability[END_REF]. Then, one may construct a function P ε satisfying y ε (•, t) -P ε (•, t) L 2 (0,1) ≤ c ε 1/2 , ∀t ∈ [0, T ],

for some positive constant c > independent of ε.

Compared to [START_REF]Asymptotic analysis of an advection-diffusion equation involving interacting boundary and internal layers[END_REF] where a rate of order O(ε 3/2 ) is obtained for the constant case, we thus get a rate of order O(ε 1/2 ). As we shall see, this is due to the fact that the non constant case under consideration here involves less explicit computations and developments leading to a more tedious analysis. The outline is as follows. In Section 2, we set up the matched asymptotic method by considering two terms in the expansion. In particular, we employ fundamental solutions for a parabolic equation with unbounded coefficients to construct an approximation of the solution in the internal layers of size ε 1/2 (we refer to Lemma 2.2). The method leads to a so-called composite approximations P ε given by (32). Then, in section 3 by considering the system satisfied by y ε -P ε and using a Gronwall-Bellman type inequality (see Lemma 3.4), we prove an a priori estimate for y ε -P ε leading to Theorem 1.1. Some computations are collected in the appendice. In Section 4, we consider some explicit simple examples of transport velocity M and evaluate numerically the norm y ε -P ε L ∞ (0,T ;L 2 (0,1)) with respect to ε in agreement with our theoretical estimate. We conclude with some remarks in Section 5.

Matched asymptotic expansions and approximate solutions

In order to construct an asymptotic approximation of the solution y ε of (1), we use the method of matched asymptotic expansions ( [START_REF] Kevorkian | Multiple scale and singular perturbation methods[END_REF][START_REF] Hubert | Vibration and coupling of continuous systems[END_REF][START_REF] Van Dyke | Perturbation methods in fluid mechanics[END_REF][START_REF] Eckhaus | Asymptotic analysis of singular perturbations[END_REF][START_REF] Ilin | Matching of asymptotic expansions of solutions of boundary value problems[END_REF]). We also refer to [START_REF] Gie | Singular perturbations and boundary layers[END_REF] for a recent exposition. The solution y ε exhibits two inner regions: an internal layer located along the characteristic {(x, t) ∈ Q T , x -X(t; 0, 0) = 0} and a boundary layer living along x = 1. The internal layer is of size O(ε 1/2 ) while the boundary layer is of size O(ε). The outer region is the subset of (0, 1) consisting of the points far from the internal and boundary layers, it is of O(1) size. The occurrence of these three distinct regions require to introduce three distinct asymptotic expansions. The first one, the so-called outer expansion, lives far away from the inner regions and is given by

m k=0 ε k y k (x, t), (x, t) ∈ Q T , x -X(t; 0, 0) = 0, x < 1,
for some m ∈ N . A second one, the so-called first inner expansion, living in the neighborhood of {(x, t) ∈ Q T , x -X(t; 0, 0) = 0} is given by

m k=0 ε k/2 W k/2 (w, t), w := x -X(t; 0, 0) ε 1/2 ∈ - X(t; 0, 0) ε 1/2 , 1 -X(t; 0, 0) ε 1/2 , t ∈ (0, T ).
Last, a third one, the so-called second inner expansion, living along x = 1, is given by

m k=0 ε k/2 Y k/2 (z, τ, t), z := 1 -x ε ∈ (0, ε -1 ), τ := T 1 -t ε 1/2 ,
where T 1 is the unique solution of (3). In particular, these expansions make appear several variables, at different scales, namely, x, t, z = ε -1 (1 -x), w = ε -1/2 (x -X(t; 0, 0)) and τ = ε -1/2 (T 1 -t).

We will construct outer and inner expansions which will be valid in the so-called outer and inner regions, respectively. There are intermediate regions between the outer region and the inner regions, with size O(ε γ ), γ ∈ (0, 1). To construct an approximate solution we require that inner and outer expansions coincide in each intermediate region, then some conditions must be satisfied in that region by the corresponding inner and outer expansions. These conditions are the so-called matching asymptotic conditions.

The strategy is as follows. We first identify the functions y k , k = 0, . . . , m, in the outer region. Then, we identify the functions W k/2 , k = 0, . . . , m, of the first inner expansion satisfying the matching conditions (with the y k ). This allows to define an expansion in the form p ε = m k=0 ε k/2 p k/2 , valid far away from x = 1, as a linear combination of the functions y k and W k/2 . Then, we identify the functions Y k , k = 0, . . . , m, of the second inner expansion satisfying the matching conditions (with the p k/2 ). Eventually, we define an expansion P ε , valid in the whole domain Q T , as a linear combination of the functions p k/2 and Y k/2 , and supposed to be an approximation of y ε . In this work, for simplicity we restrict ourselves to the case m = 1 as it will allow to get an approximation with rate ε 1/2 .

A property of the characteristics

We assume that M is a positive and smooth function in Q T . Let t → X(t; x, s) be the characteristic through (x, s) ∈ Q T defined as the solution in Q T of equation [START_REF]On the controllability of an advection-diffusion equation with respect to the diffusion parameter: asymptotic analysis and numerical simulations[END_REF]. Using a classical result of differentiation with respect to a parameter of solutions of a differential equation we get that X is a smooth function of the variables t, x, and s. In particular,

∂ t X x = M x (X, t)X x , X x (s; x, s) = 1,
i.e. X x satisfies the integral equation

X x (t; x, s) = 1 + t s M x (X(σ; x, s), σ)X x (σ; x, s) dσ.
We set a(t) := X(t; 0, 0), t ≥ 0.

Let γ 0 (x, t) denote the unique solution of the equation X(γ 0 (x, t); x, t) = 0, for x < a(t).

Differentiating equality (4) with respect to x gives X t (γ 0 (x, t); x, t)γ 0x (x, t) + X x (γ 0 (x, t); x, t) = 0 for all (x, t) with x < a(t), hence

γ 0x (x, t) = - X x (γ 0 (x, t); x, t) M (0, γ 0 (x, t)) , for x < a(t), (5)
since X t (γ 0 (x, t); x, t) = M (X(γ 0 (x, t); x, t), γ 0 (x, t)) = M (0, γ 0 (x, t)).

Outer expansion

Putting y 0 (x, t) + εy 1 (x, t) into equation (1) 1 , the identification of the powers of ε yields

ε 0 : y 0 t + M y 0 x = 0, ε : y 1 t + M y 1 x = y 0 xx . (6)
Taking the initial and boundary conditions into account we define y 0 and y 1 as functions satisfying the transport equations, respectively,

     y 0 t (x, t) + M (x, t)y 0 x (x, t) = 0, (x, t) ∈ Q T , y 0 (0, t) = v(t), t ∈ (0, T ), y 0 (x, 0) = y 0 (x), x ∈ (0, 1),      y 1 t (x, t) + M (x, t)y 1 x (x, t) = y 0 xx (x, t), (x, t) ∈ Q T , y 1 (0, t) = 0, t ∈ (0, T ), y 1 (x, 0) = 0, x ∈ (0, 1).
We find explicit representations of solutions of these equations by using the method of characteristics. We get

y 0 (x, t) = y 0 (X(0; x, t), x > a(t), v (γ 0 (x, t)) , x < a(t), (7)
and

y 1 (x, t) =          t 0 y 0 xx (X(s; x, t), s)ds, x > a(t), t γ0(x,t)
y 0 xx (X(s; x, t), s)ds, x < a(t).

Inner expansion along the characteristic x -a(t) = 0

We consider the change of variable w = x-a(t) ε 1/2 and function W ε (w, t) = y ε (x, t). W ε satisfies the equation

W ε t (w, t) + M (ε 1/2 w + a(t), t) -M (a(t), t) ε 1/2 W ε w (w, t) -W ε ww (w, t) = 0. ( 8
)
Using the Taylor expansion

M ε 1/2 w + a(t), t -M (a(t), t) = ε 1/2 wM x (a(t), t) + εw 2 2 M xx (a(t), t) + O(ε 3/2 ),
then, putting W 0 (w, t) + εW 1/2 (w, t) into equation ( 8), the identification of the powers of ε yields

ε 0 : W 0 t (w, t) + M x (a(t), t)w W 0 w (w, t) -W 0 ww (w, t) = 0, ε 1/2 : W 1/2 t (w, t) + M x (a(t), t)w W 1/2 w (w, t) -W 1/2 ww (w, t) = - w 2 2 M xx (a(t), t)W 0 w (w, t).
Obviously, the main difference with respect to the case M constant considered in [START_REF]Asymptotic analysis of an advection-diffusion equation involving interacting boundary and internal layers[END_REF] is the occurence of some aditionnal unbounded terms with respect to the variable w.

To get the asymptotic matching conditions we write that, for any fixed t and large w,

W 0 (w, t) + ε 1/2 W 1/2 (w, t) = y 0 (x, t) + εy 1 (x, t) + O(ε 2 ).
Rewriting the right-hand side of the above equality in terms of w, t, and using Taylor expansions we have

W 0 (w, t) + ε 1/2 W 1/2 (w, t) + εW 1 (w, t) = y 0 ε 1/2 w + a(t), t + εy 1 (ε 1/2 w + a(t), t) + O(ε 2 ) = y 0 (a(t), t) + ε 1/2 wy 0 x (a(t), t) + εw 2 2 y 0 xx (a(t), t) + O(ε 3/2 ).
Therefore, the matching conditions read

W 0 (w, t) ∼ y 0 ((a(t)) ± , t), W 1/2 (w, t) ∼ y 0 x ((a(t)) ± , t)w, as w → ±∞. (9)
Consequently, we obtain that the function

W 0 must satisfy    W 0 t (w, t) + M x (a(t), t)w W 0 w (w, t) -W 0 ww (w, t) = 0, (w, t) ∈ R × (0, T ), lim w→±∞ W 0 (w, t) = lim x→(a(t)) ± y 0 (x, t), t ∈ (0, T ). ( 10
)
In view of (7), we have lim x→(a(t)) + y 0 (x, t) = y 0 (0) and lim x→(a(t)) -y 0 (x, t) = v(0). Remark that these limits do not depend on t. Similarly, the function W 1/2 must satisfy

       W 1/2 t (w, t) + M x (a(t), t)w W 1/2 w (w, t) -W 1/2 ww (w, t) = - w 2 2 M xx (a(t), t)W 0 w (w, t), (w, t) ∈ R × (0, T ), lim w→±∞ W 1/2 (w, t) -y 0 x (a(t) ± , t) w = 0, t ∈ (0, T ).
(11) In view of (7), y 0

x (a(t) + , t) = (y 0 ) (1) (0) X x (0; a(t), t) and y 0 x (a(t) -, t) = v (1) (0) γ 0x (a(t) -, t) where γ 0x is given by [START_REF] Aronson | Parabolic equations with unbounded coefficients[END_REF]. Remark that these limits do depend on the variable t.

2.4 Representation of the functions W 0 and W 1/2 We express the functions W 0 and W 1/2 in term of the fundamental solution of the heat equation.

Fundamental solutions

Let us consider the differential operator

LU (w, t) := U t (w, t) + M x (a(t), t)w U w (w, t) -U ww (w, t) in R × (0, T ). ( 12
)
Equation LU = 0 is a parabolic equation with an unbounded coefficient. Linear second order parabolic equations with smooth but unbounded coefficients are studied in [START_REF] Aronson | Parabolic equations with unbounded coefficients[END_REF] (see also [START_REF] Ilin | Second-order linear equations of parabolic type[END_REF][START_REF] Friedman | Partial differential equations of parabolic type[END_REF][START_REF] Bodanko | Sur le problème de Cauchy et les problèmes de Fourier pour les équations paraboliques dans un domaine non borné[END_REF]). In [START_REF] Aronson | Parabolic equations with unbounded coefficients[END_REF], the authors consider a second order linear differential operator in the form

Lu(w, t) = u t (w, t) -c ij (w, t)u wiwj (w, t) -c i (w, t)u wi (w, t) -c(w, t)u(w, t), (w, t) ∈ R n × (0, T ),
and prove the existence of a so-called fundamental solution for Lu = 0. Our operator L is obviously of the form of the operator L with n = 1, c ≡ 0, c 11 ≡ 0 and c 1 (w, t) = M x (a(t), t)w.

Definition 2.1. A function K(w, t; ξ, s) defined for w, ξ ∈ R n and 0 ≤ s < t ≤ T is said to be a fundamental solution of Lu = 0 if it has the following properties:

• considered as a function of (w, t) for each fixed (ξ, s) ∈ R n × [0, T ], the derivatives of K which appear in L exist and are continuous;

• L K = 0 in R n × (s, T ]; • If g = g(w)
is a continuous function with compact support in R n then lim (w,t)→(w0,s + ) R n K(w, t; ξ, s)g(ξ) dξ = g(w 0 ).

We shall employ the following lemma.

Lemma 2.1. Let f = f (w, t) be a given function in R n × (0, T ). Assume that f is Hölder continuous on every compact subset of R n × (0, T ). Then

u(w, t) = R n K(w, t; ξ, 0)g(ξ) dξ + t 0 R n K(w, t; ξ, τ )f (ξ, τ ) dξ dτ is a solution of the Cauchy problem Lu(w, t) = f (w, t) in R n × (0, T ), u(w, 0) = g(w) in R n .
In the particular case of the operator L, the solution of the corresponding Cauchy problem has an explicit representation in term of the fundamental solution of the heat equation

H s -H vv = 0 in R × (0, T ) defined as follows H(v, s; ξ, τ ) = H 0 (v -ξ; s -τ ), v, ξ ∈ R, 0 ≤ τ < s ≤ T, (13)
with

H 0 (v, s) := 1 √ 4πs e -v 2 4s , v ∈ R, s ∈ (0, T ). ( 14
)
Precisely, we have the following result.

Lemma 2.2. Let L be the parabolic operator defined by [START_REF] Friedman | Partial differential equations of parabolic type[END_REF] and H be the fundamental solution of the heat equation defined by [START_REF] Gie | Singular perturbations and boundary layers[END_REF]. Then: i) Equation LU = 0 has a fundamental solution given by

K(w, t; ξ, τ ) = 1 β(τ ) H w β(t) , B(t); ξ β(τ ) , B(τ ) , w, ξ ∈ R, 0 ≤ τ < t ≤ T, (15)
with α(t) := M x (a(t), t), β(t) := e t 0 α(s) ds , B(t

) := t 0 1 β(s) 2 ds.
ii) Let f = f (w, t) be a continuous function in R×[0, T ] and Hölder continuous with respect to w uniformly for t ∈ [0, T ]. Let g = g(w) be a continuous function with compact support in R. Then,

U (w, t) = R K(w, t; ξ, 0)g(ξ) dξ + t 0 R K(w, t; ξ, τ )f (ξ, τ ) dξdτ ( 16
)
is a solution of the Cauchy problem

LU (w, t) = f (w, t) in R × (0, T ), U (w, 0) = g(w) in R. (17)
Proof. We have

β (t) = α(t)β(t), B (t) = 1 β(t) 2 , β(0) = 1, B(0) = 0.
Note also that B is a stricty increasing function so that 0 ≤ τ < t ≤ T is equivalent to 0 ≤ B(τ ) < B(t) ≤ B(T ). Consider the function K defined by [START_REF] Ilin | Matching of asymptotic expansions of solutions of boundary value problems[END_REF]. We have

K(w, t; ξ, τ ) = 1 β(τ ) 4π(B(t) -B(τ )) e - ( w β(t) - ξ β(τ ) ) 2 4(B(t)-B(τ )) , 0 ≤ τ < t ≤ T.
For ξ, τ fixed, ξ ∈ R, 0 ≤ τ < t ≤ T , we have by a direct calculation

LK(w, t; ξ, τ ) = 1 β(τ )β 2 (t) H 0s w β(t) - ξ β(τ ) , B(t) -B(τ ) -H 0vv w β(t) - ξ β(τ ) , B(t) -B(τ ) .
Then, LK(w, t; ξ, τ ) = 0 for any w ∈ R and t ∈ (0, T ). If g = g(w) is a continuous function with compact support in R we have R K(w, t; ξ, τ )g(ξ) dξ = 1

β(τ ) 4π(B(t) -B(τ )) R e - ( w β(t) - ξ β(τ ) ) 2 4(B(t)-B(τ )) g(ξ) dξ.
Using the change of variable

ξ β(τ ) -w β(t) 2 √ (B(t)-B(τ )) = v, it holds that R K(w, t; ξ, τ )g(ξ) dξ = 1 √ π R e -v 2 g β(τ ) β(t) w + 2β(τ ) B(t) -B(τ ) v dv,
and deduce that lim t→τ + R K(w, t; ξ, τ )g(ξ) dξ = g(w) for any w ∈ R.

According to Definition 2.1, it follows that the function K given by ( 15) is a fundamental solution of the equation LU = 0. Moreover, using that β(0) = 1 and B(0) = 0, we have

R K(w, t; ξ, 0)g(ξ) dξ = 1 √ π R e -v 2 g w β(t) + 2 B(t) v dv,
implying that lim t→0 + R K(w, t; ξ, 0)g(ξ) dξ = g(w).

We also have, for any ξ fixed, LK(w, t; ξ, 0) = 0 from which we check that L R K(w, t; ξ, 0)g(ξ) dξ = 0. Consider now the second term in the right-hand side of ( 16) which we denote Ũ (w, t). We have

L( Ũ )(w, t) = lim τ →t - R K(w, t; ξ, τ )f (ξ, τ ) dξ + t 0 R L(K)(w, t; ξ, τ )f (ξ, τ ) dξdτ = lim τ →t - R K(w, t; ξ, τ )f (ξ, τ ) dξ.
Using again the change of variable

ξ β(τ ) -w β(t) 2 √ (B(t)-B(τ )) = v, it holds that lim τ →t - R K(w, t; ξ, τ )f (ξ, τ ) dξ = 1 √ π lim τ →t - R e -v 2 f β(τ ) β(t) w + 2β(τ ) (B(t) -B(τ )) v, τ dv = f (w, t).
We conclude that U is a solution of the Cauchy problem [START_REF] Kevorkian | Multiple scale and singular perturbation methods[END_REF].

We now introduce the functions g k/2 = g k/2 (w, t), k = 0, 1, defined as functions on R × (0, T ) and satisfying lim w→±∞ g 0 (w, t) -y 0 (a(t) ± , t) = 0, lim

w→±∞ g 1/2 (w, t) -y 0 x (a(t) ± , t) w = 0, t ∈ (0, T ).
We also introduce the functions

f k/2 = f k/2 (w, t), k = 0, 1, defined on R × (0, T ) by f 0 (w, t) := 0, f 1/2 (w, t) := - w 2 2 M xx (a(t), t) W 0 w (w, t). (18)

Explicit expression of W 0

In order to apply the second part of Lemma 2.2, it remains to choose an initial condition g (for the boundary value problem satisfied by W 0 ) so that the corresponding solution W 0 given generically by [START_REF] Ilin | Second-order linear equations of parabolic type[END_REF] satisfies the asymptotic conditions as w → ±∞ (see [START_REF] Deuring | L 2 -stability independent of diffusion for a finite elementfinite volume discretization of a linear convection-diffusion equation[END_REF]), i.e. lim w→∞ W 0 (w, t) = y 0 (0) and lim w→-∞ W 0 (w, t) = v(0). The simplest choice, discussed in [START_REF]Asymptotic analysis of an advection-diffusion equation involving interacting boundary and internal layers[END_REF], is given by g = g 0 with g 0 (w) := y 0 (0), w ≥ 0,

v(0), w < 0, (19)
leading to the explicit expression:

W 0 (w, t) = y 0 (0) -v(0) 2 erf w 2β(t) B(t) + y 0 (0) + v(0) 2 . ( 20
)
Recall that erf is the error function defined by erf (s) = 2 √ π s 0 e -v 2 dv for all s ∈ R, and that the complementary error function is defined by erf c(s

) = 1 -erf (s) = 2 √ π
+∞ s e -v 2 dv. Using the asymptotic behavior of the error function, we check that W 0 satisfies the prescribed asymptotic behavior in [START_REF] Deuring | L 2 -stability independent of diffusion for a finite elementfinite volume discretization of a linear convection-diffusion equation[END_REF]. Remark also that lim t→0 + W 0 (w, t) equals y 0 (0) if w ≥ 0 and v 0 (0) if w < 0.

Integral representation of W 1/2

The computation of W 1/2 solution of the boundary value problem [START_REF] Eckhaus | Asymptotic analysis of singular perturbations[END_REF] posed over R × (0, T ) is more delicate since, first there is a right side, and second, because the asymptotic limits as w → ∞ depend on the time variable, while the function g in ( 17) not. This fact requires a special treatment.

Let us consider the function g 1/2 (w, t) = y 0 x (a(t) ± , t) w. We have, for (w,

t) ∈ (] -∞, 0[∪]0, +∞[)×(0, T ), L(g 1/2 )(w, t) = w y 0 xt (a(t) ± , t) + a (t)y 0 xx (a(t) ± , t) + M x (a(t), t)wy 0 x (a(t) ± , t) = w y 0 xt (a(t) ± , t) + M (a(t), t)y 0 xx (a(t) ± , t) + M x (a(t), t)wy 0 x (a(t) ± , t).
Differentiating equation ( 6) 1 with respect to x yields

y 0 xt (x, t) + M x (x, t)y 0 x (x, t) + M (x, t)y 0 xx (x, t) = 0, (x, t) ∈ Ω + ∪ Ω -, with Ω + := {(x, t) ∈ Q T , x > a(t)} and Ω -:= {(x, t) ∈ Q T , x < a(t)}. Letting x → a(t) ± , we get y 0 xt (a(t) ± , t) + M x (a(t), t)y 0 x (a(t) ± , t) + M (a(t), t)y 0 xx (a(t) ± , t) = 0. Therefore, for (w, t) ∈ (] -∞, 0[∪]0, +∞[) × (0, T ), we have L(g 1/2 )(w, t) = -wM x (a(t), t)y 0 x (a(t) ± , t) + wM x (a(t), t)y 0 x (a(t) ± , t) = 0.
When y 0 x (a(t) ± , t) has a jump on w = 0, the function (g 1/2 ) ww becomes singular on w = 0. This leads to regularize the function

g 1/2 as follows. Let ρ ∈ C ∞ (R) such that ρ(w) = 1 for |w| ≥ 2, 0 for |w| ≤ 1.
We define

g 1/2 (w, t) := ρ(w) g 1/2 (w, t), w ∈ R, t ∈ (0, T ). ( 21
)
Clearly, g 1/2 ∈ C ∞ (R), g 1/2 (w, t) = g 1/2 (w, t) for |w| ≥ 2 and g 1/2 (w, t) = 0 for |w| ≤ 1. Define

g 1/2 (w, t) = L( g 1/2 )(w, t), w ∈ R, t ∈ (0, T ), (22)
and then W 1/2 as follows

W 1/2 (w, t) = g 1/2 (w, t) + t 0 R K(w, t; ξ, τ ) f 1/2 (ξ, τ ) -g 1/2 (ξ, τ ) dξdτ, w ∈ R, t ∈ (0, T ). ( 23
) Lemma 2.3. The function W 1/2 defined by (23) satisfies (11), namely L(W 1/2 )(w, t) = f 1/2 (w, t) in R×(0, T ) and lim w→±∞ W 1/2 (w, t) -g 1/2 (w, t) = 0 for all t in (0, T ). Proof. According to Lemma 2.2 we have L(W 1/2 )(w, t) = f 1/2 (w, t) in R × (0, T ).
It remains to check the asymptotic condition as w → ±∞. For any (w, t) ∈ R × (0, T ), let us consider the integral

I 1 (w, t) := t 0 R K(w, t; ξ, τ )f 1/2 (ξ, τ ) dξdτ.
In view of ( 18) and ( 20), we obtain explicitly

f 1/2 (w, t) = - w 2 2 M xx (a(t), t) y 0 (0) -v(0) 2 2 √ π e - w 2 4β(t) 2 B(t) 1 2β(t) B(t) , then I 1 (w, t) = y 0 (0) -v(0) 2 t 0 R 1 β(τ ) 4π(B(t) -B(τ )) e - ( w β(t) - ξ β(τ ) ) 2 4(B(t)-B(τ )) × × - ξ 2 2 M xx (a(τ ), τ ) 2 √ π e - ξ 2 4β(τ ) 2 B(τ ) 1 2β(τ ) B(τ ) dξdτ.
For the convergence of the integral I 1 (ω, t) in the neighborhood of τ = 0, we consider the integral I 1 (ω, t ), with 0 < t < t fixed, and use the change of variable

ξ 2β(τ ) √ B(τ ) = v. This gives I 1 (w, t ) := y 0 (0) -v(0) 2 t 0 R 1 β(τ ) 4π(B(t) -B(τ )) e - ( w β(t) -2 √ B(τ )v ) 2 4(B(t)-B(τ )) × × -2β(τ ) 2 B(τ )v 2 M xx (a(τ ), τ ) 2 √ π e -v 2 dvdτ.
The integrand is clearly uniformly bounded in w by c(t , t)e -v 2 with a constant c(t , t) depending only on t and t. Moreover, for any fixed v and τ , the integrand tends to 0 as |w| → +∞. Using a result of limit passage on integrals depending on a parameter we get lim |w|→+∞ I 1 (w, t ) = 0.

For the convergence of the integral I 1 (w, t) in the neighborhood of τ = t, we consider the integral

I 2 (w, t) := y 0 (0) -v(0) 2 t t R 1 β(τ ) 4π(B(t) -B(τ )) e - ( w β(t) - ξ β(τ ) ) 2 4(B(t)-B(τ )) × × - ξ 2 2 M xx (a(τ ), τ ) 2 √ π e - ξ 2 4β(τ ) 2 B(τ ) 1 2β(τ ) B(τ ) dξdτ.
Using the change of variable

ξ β(τ ) -w β(t) 2 √ B(t)-B(τ )
= v, we get

I 2 (w, t) = y 0 (0) -v(0) 2 1 √ π t t R e -v 2   - β(τ ) β(t) w + 2β(τ ) B(t) -B(τ ) v 2 2    M xx (a(τ ), τ )× × 2 √ π e - ( β(τ ) β(t) w+2β(τ ) √ B(t)-B(τ ) v ) 2 4β(τ ) 2 B(τ ) 1 2β(τ ) B(τ ) dvdτ
and conclude as before that lim |w|→+∞ I 2 (w, t) = 0 and then that lim |w|→+∞ I 1 (w, t) = 0. Consider now the integral

I 3 (w, t) : = t 0 R K(w, t; ξ, τ )g 1/2 (ξ, τ ) dξdτ = t 0 R 1 β(τ ) 4π(B(t) -B(τ )) e - ( w β(t) - ξ β(τ ) ) 2 4(B(t)-B(τ )) g 1/2 (ξ, τ ) dξdτ.
Clearly, g 1/2 is a smooth function on R × (0, T ) and satisfies

g 1/2 (w, t) = 0 for |w| ≥ 2. The change of variable ξ β(τ ) -w β(t) 2 √ (B(t)-B(τ )) = v gives I 3 (w, t) = 1 √ π t 0 R e -v 2 g 1/2 β(τ ) β(t) w + 2β(τ ) B(t) -B(τ ) v, τ dvdτ.
The integrand is uniformly bounded in w by c(t)e -v 2 , where c(t) depends only on t. Moreover, the integrand tends to 0 as |w| → +∞. Then lim |w|→+∞ I 3 (w, t) = 0.

We then have

lim |w|→+∞ t 0 R K(w, t; ξ, τ ) f 1/2 (ξ, τ ) -g 1/2 (ξ, τ ) dξdτ = 0.
It results from ( 23) that, for each fixed t, lim w→±∞ W 1/2 (w, t)-g 1/2 (w, t) = 0. Since lim w→±∞ g 1/2 (w, t)g 1/2 (w, t) = 0, we conclude that lim w→±∞ W 1/2 (w, t) -g 1/2 (w, t) = 0.

Composite asymptotic approximation outside the boundary layer

We now define a composite approximation of the solution far from the x = 1 obtained by adding, at each order, the inner and outer expansions and then by subtracting their common part. At the first order, the common part of y 0 (x, t) and W 0 (w, t) (defined by ( 7) and ( 20) respectively) is equal to y 0 (0) for x > a(t) and to v(0) for x < a(t). Thus, the first term of the composite approximation outside the boundary layer along x = 1 is given by

p 0 (x, t) = y 0 (x, t) + W 0 (w, t) -y 0 (a(t) ± , t). ( 24
)
We check that the function p 0 is continuous along the characteristic as it satisfies lim x-a(t)→0 ± p 0 (x, t) = W 0 (0, t). The second term of the composite approximation is given by

p 1/2 (x, t) = W 1/2 (w, t) -y 0 x (a(t) ± , t), (25)
where W 1/2 is defined by [START_REF] Shih | A novel uniform expansion for a singularly perturbed parabolic problem with corner singularity[END_REF]. Clearly, p 1/2 is also continuous along the characteristic. Then the following quantity is defined to be an asymptotic approximation of y ε , outside the boundary layer along x = 1,

p ε (x, t) := p 0 (x, t) + ε 1/2 p 1/2 (x, t), (x, t) ∈ Q T .
We easily verify the following property.

Proposition 2.1. Assume that v ∈ C 1 ([0, T ]) and y 0 ∈ C 1 ([0, 1]). Then p ε belongs to C 1 ([0, 1] × (0, T ]).

Inner expansion along x = 1

We consider the change of variables z = 1-x ε , τ = T1-t ε 1/2 , and function Y ε (z, τ, t) = y ε (x, t). Recall that T 1 defined by a(T 1 ) = 1 is the time at which the characteristic x -a(t) = 0 starting from the point (0, 0) intersects the right extremity of the domain. We assume that T 1 exists. Moreover, since M > 0, T 1 is unique. The function Y ε satisfies the equation

Y ε t (z, τ, t) - 1 ε 1/2 Y ε τ (z, τ, t) - M (1 -εz, t) ε Y ε z (z, τ, t) - 1 ε Y ε zz (z, τ, t) = 0. ( 26
)
Using the Taylor expansion 26), the identification of the powers of ε yields

M (1 -εz, t) = M (1, t) -εzM x (1, t) + ε 2 z 2 2 M xx (1, t) - ε 3 z 3 3! M xxx (1, t) + • • • , and putting Y 0 (z, τ, t) + ε 1/2 Y 1/2 (z, τ, t) into equation (
ε -1 : Y 0 zz (z, τ, t) + M (1, t)Y 0 z (z, τ, t) = 0, ε -1/2 : Y 1/2 zz (z, τ, t) + M (1, t)Y 1/2 z (z, τ, t) = -Y 0 τ (z, τ, t).
It is important to note that the functions Y k depend on three variables, namely z, t but also τ . As it is standard, the variable z = (1 -x)/ε is introduced to describe the boundary layer at x = 1 -. Here, the variable τ = T1-t ε 1/2 allows to take into account the interaction of the internal and boundary layers by making a zoom around the point (1, T 1 ). If we do not introduce this variable τ , we see notably that Y 1/2 solves the same ordinary differential equation than Y 0 , and the analysis leads to an error estimate (for the L ∞ (L 2 ) norm) of order ε 1/4 only. Remark that the variable T 1 , as a solution of a(T 1 ) = 1 is not explicit. Therefore, instead of τ = T1-t ε 1/2 , we could have introduce τ = 1-a(t) ε 1/2 ; however, we have observed that this leads to some incompatibilities with the matching condition.

We impose Y k/2 (0, τ, t) = 0 for k = 0, 1. To get the asymptotic matching conditions we write that, for any fixed τ, t and large z,

Y 0 (z, τ, t) + ε 1/2 Y 1/2 (z, τ, t) = p 0 (x, t) + ε 1/2 p 1/2 (x, t) + O(ε 2 ).
In order to identify at each order the appropriate matching conditions, we need to rewrite the right-hand side of the above equality in terms of z and τ , t being fixed. We have the following equalities

x = 1 -εz, w = x-a(t) ε 1/2 = -ε 1/2 z + 1-a(t) ε 1/2 , a(T 1 ) = 1, a (t) = X t (t; 0, 0) = M (a(t), t), a (T 1 ) = M (1, T 1 ). Writing a(t) = a T 1 -ε 1/2 τ , M (1, T 1 ) = M 1, t + ε 1/2 τ
, and using Taylor expansions we have

a(t) = a(T 1 ) -ε 1/2 τ a (T 1 ) + O(ε) = 1 -ε 1/2 τ M (1, T 1 ) + O(ε), M (1, T 1 ) = M (1, t) + O(ε 1/2 ), then 1 -a(t) ε 1/2 = τ M (1, t) + O(ε 1/2 ). It results that w = -ε 1/2 z + τ M (1, t) + O(ε 1/2 ). Writing y 0 (x, t) = y 0 (1 -εz, t), W 0 (w, t) = W 0 -ε 1/2 z + τ M (1, t) + O(ε 1/2 ), t , and W 1/2 (w, t) = W 1/2 -ε 1/2 z + τ M (1, t) + O(ε 1/2
), t , and using again Taylor expansions we have

p 0 (x, t) = y 0 (x, t) + W 0 (w, t) -y 0 (a(t) ± , t) = y 0 (1, t) + W 0 (τ M (1, t), t) -y 0 (a(t) ± , t) -ε 1/2 zW 0 w (τ M (1, t), t) + O(ε 1/2 ),
and

ε 1/2 p 1/2 (x, t) =ε 1/2 W 1/2 (w, t) -y 0 x (a(t) ± , t)w =ε 1/2 W 1/2 (τ M (1, t), t) -τ M (1, t) y 0 x (a(t) ± , t) + O(ε).
We deduce that p 0 (x, t)

+ ε 1/2 p 1/2 (x, t) = C 0 (τ, t) + ε 1/2 C 1/2 (z, τ, t) + O(ε), with C 0 (τ, t) =y 0 (1, t) + W 0 (τ M (1, t), t) -y 0 (a(t) ± , t), C 1/2 (z, τ, t) =W 1/2 (τ M (1, t), t) -τ M (1, t)y 0 x (a(t) ± , t) -zW 0 w (τ M (1, t), t) . (27)
Therefore, the matching conditions read

Y 0 (z, τ, t) ∼ C 0 (τ, t), Y 1/2 (z, τ, t) ∼ C 1/2 (z, τ, t), as z → ∞. (28)
• We define Y 0 as a solution of

   Y 0 zz (z, τ, t) + M (1, t)Y 0 z (z, τ, t) = 0, (z, τ, t) ∈ R + × R × (0, T ), Y 0 (0, τ, t) = 0, lim z→+∞ Y 0 (z, τ, t) = C 0 (τ, t), t ∈ (0, T ).
The solution is

Y 0 (z, τ, t) = C 0 (τ, t) 1 -e -M (1,t)z , (z, τ, t) ∈ R + × R × [0, T ].
• We define Y 1/2 as a solution of

   Y 1/2 zz (z, τ, t) + M (1, t)Y 1/2 z (z, τ, t) = -Y 0 τ (z, τ, t), (z, τ, t) ∈ R + × R × (0, T ), Y 1/2 (0, τ, t) = 0, lim z→+∞ Y 1/2 (z, τ, t) -C 1/2 (z, τ, t) = 0, t ∈ (0, T ). Writing that Y 0 τ (z, τ, t) = C 0,τ (τ, t) 1 -e -M (1,t)z , we obtain, for (z, τ, t) ∈ R × R × (0, T ), Y 1/2 (z, τ, t) = C 1/2 (z, τ, t) + e -M (1,t)z -C 1/2 (0, τ, t) -W 0 w (τ M (1, t), t)z , (z, τ, t) ∈ R + × R × [0, T ].
(29)

Asymptotic composite approximation in Q T

We are now in position to define what is supposed to be an asymptotic approximation of the solution y ε . We proceed as before by adding at each order the function p k/2 , approximation outside the boundary layer along x = 1, and the function Y k/2 , approximation in the boundary layer along x = 1, then subtracting their common part. At the first order, the composite approximation is given by

P 0 (x, t) = p 0 (x, t) + Y 0 (z, τ, t) -C 0 (τ, t) = p 0 (x, t) -C 0 (τ, t)e -M (1,t)z , (30)
i.e. explicitly

P 0 (x, t) = y 0 (x, t) + W 0 (w, t) -y 0 (a(t) ± , t) -y 0 (1, t) + W 0 (τ M (1, t), t) -y 0 (a(t) ± , t) e -M (1,t)z .
Repeating the arguments at the next order, we define

P 1/2 (x, t) = p 1/2 (x, t) + Y 1/2 (z, τ, t) -C 1/2 (z, τ, t), (31)
then we define an asymptotic composite approximation of y ε in Q T by

P ε (x, t) := P 0 (x, t) + ε 1/2 P 1/2 (x, t), (x, t) ∈ Q T . ( 32
)
3 Convergence of the sequence (P ε ) (ε>0) -Proof of Theorem 1.1

This section is devoted to the study of the convergence of the sequence (P ε ) ε>0 stated in Theorem (1.1).

In order to prove this theorem we need to establish a number of preliminary results. We define the error as follows:

z ε (x, t) := P ε (x, t) -y ε (x, t) -θ ε (x, t), (x, t) ∈ Q T , (33)
where θ ε is the initial layer corrector defined as a solution of the equation

     θ ε t (x, t) + M (x, t)θ ε x -εθ ε xx (x, t) = 0, (x, t) ∈ Q T , θ ε (0, t) = θ ε (1, t) = 0, t ∈ (0, T ), θ ε (x, 0) = P ε (x, 0) -y ε (x, 0),
x ∈ (0, 1).

(34)

Preliminary results

From now on, in order to shorten some equations, we shall use the following notation:

L ε y := y t + M (x, t)y x -εy xx . (35)
3.1.1 Estimate of the initial layer corrector θ ε

The following lemma gives an exponential decay property of the initial layer corrector.

Lemma 3.1. Let θ ε be the solution of problem (34) and γ ∈ (0, 1/2]. There exists a constant c independent of ε such that

θ ε (•, t) L 2 (0,1) ≤ c e -M (1,0)ε γ ε + c ε 1/2 e -M 2 0 2ε γ t , ∀t ∈ [0, T ]. ( 36
)
Proof. i) We first check that the initial data θ ε (•, 0) is given by

θ ε (x, 0) = -y 0 (1)e -M (1,0)z , ∀x ∈ (0, 1], (37)
showing that the initial condition gets concentrated in the neighborhood of x = 1. Indeed, from (30)-(32) we have

θ ε (x, 0) = 1 k=0 ε k/2 lim t→0 P k/2 (x, t) -y 0 (x) = lim t→0 p 0 (x, t) -C 0 (τ, t)e -M z -y 0 (x) + lim t→0 ε 1/2 p 1/2 (x, t) + Y 1/2 (z, τ, t) -C 1/2 (z, τ, t) , x ∈ (0, 1).
We have from (30)

lim t→0 p 0 (x, t) -C 0 (τ, t)e -M (1,t)z -y 0 (x) = lim t→0 W 0 (w, t) -y 0 (0) -C 0 (τ, t)e -M (1,t)z = -lim t→0 C 0 (τ, t)e -M (1,t)z = -lim t→0 y 0 (1, t) + W 0 (τ M (1, t), t) -y 0 (a(t)) ± , t) e -M (1,t)z = -y 0 (1)e -M (1,0)z .
Using the matching conditions of W 1/2 with y 0 , we check that lim t→0 p 1/2 (x, t) = 0, for x ∈ (0, 1] and that lim t→0

Y 1/2 (z, τ, t) -C 1/2 (z, τ, t) = lim t→0 -W 1/2 (τ M (1, t), t) + τ y 0 x (a(t) ± , t) -zW 0 w (τ M (1, t), t) e -M
(1,t)z = 0, for x > 0.

It results that P 1/2 (x, 0) = 0 for all x ∈ (0, 1]. These computations lead to (37).

ii-a) We now introduce a C ∞ cut-off function X : R → [0, 1] such that X (s) = 0 if s ≤ 1 and X (s) = 1 if s ≥ 2 and define, for γ ∈ (0, 1/2], the function

X ε : [0, 1] → [0, 1] by X ε (x) = X 1-x ε γ .
The solution θ ε of the linear system (34) can be decomposed as

θ ε = θ ε,1 + θ ε,2 with      L ε θ ε,1 = 0, (x, t) ∈ Q T , θ ε,1 (0, t) = θ ε,1 (1, t) = 0, t ∈ (0, T ), θ ε,1 (x, 0) = X ε (x)θ ε (x, 0), x ∈ (0, 1), (38)
     L ε θ ε,2 = 0, (x, t) ∈ Q T , θ ε,2 (0, t) = θ ε,2 (1, t) = 0, t ∈ (0, T ), θ ε,2 (x, 0) = (1 -X ε (x))θ ε (x, 0), x ∈ (0, 1). ( 39
)
In view of the definition of X ε , we see that θ ε,1 (x, 0) = 0 for all x ≥ 1 -ε γ . Then, in view of (37), we check that there exists a constant c 1 > 0 independent of ε such that |θ ε,1 (x, 0)| ≤ c 1 e -M (1,0)ε γ ε for all x ∈ (0, 1). By a maximum principle, it follows that

|θ ε,1 (x, t)| ≤ c 1 e -M (1,0)ε γ ε , ∀(x, t) ∈ Q T . (40)
ii-b) Concerning θ ε,2 (•, 0), we check that θ ε,2 (x, 0) = 0 for all x ≤ 1-2ε γ and that θ ε,2 (•, 0

) L 2 (0,1) ≤ c 2 ε 1 2
for some constant c 2 > 0, independent of ε. We now estimate θ ε,2 (•, t) L 2 (0,1) . For that, for any α > 0, we check that the function ρ ε (x, t) := e -M 0 αx 2ε θ ε,2 (x, t) with M 0 := min

(x,t)∈Q T M (x, t) > 0 solves          ρ ε t -ερ ε xx + (M -αM 0 ) ρ ε x - M 0 4ε α 2 M 0 -2αM ρ ε = 0, (x, t) ∈ Q T , ρ ε (0, t) = ρ ε (1, t) = 0, t ∈ (0, T ), ρ ε (x, 0) = e -M 0 αx 2ε θ ε,2 0 (x), x ∈ (0, 1), (41)
and satisfies the estimates

d dt ρ ε (•, t) 2 L 2 (0,1) + 2ε ρ ε x (•, t) 2 L 2 (0,1) = 1 0 M x + M 0 2ε (α 2 M 0 -2αM (ρ ε (•, t)) 2 dx ≤ 1 0 C + M 2 0 2ε (α 2 -2α (ρ ε (•, t)) 2 dx with C = max (x,t)∈Q T |M x (x, t)|. Gronwall inequality leads to ρ ε (•, t) L 2 (0,1) ≤ ρ ε (•, 0) L 2 (0,1) e C 2 t e M 2 0 4ε (α 2 -2α)t ,
which is equivalent to

e -M 0 αx 2ε θ ε,2 (•, t) L 2 (0,1) ≤ e -M 0 αx 2ε θ ε,2 (•, 0) L 2 (0,1) e C 2 t e M 2 0 4ε (α 2 -2α)t .
Consequently,

θ ε,2 (•, t) L 2 (0,1) = e M 0 αx 2ε e -M 0 αx 2ε θ ε,2 (•, t) L 2 (0,1) ≤ e M 0 αx 2ε L ∞ (0,1) e -M 0 αx 2ε θ ε,2 (•, t) L 2 (0,1) ≤ e M 0 αx 2ε L ∞ (0,1) e -M 0 αx 2ε θ ε,2 (•, 0) L 2 (0,1) e C 2 t e M 2 0 4ε (α 2 -2α)t ≤ e M 0 αx 2ε L ∞ (0,1) e -M 0 αx 2ε θ ε,2 (•, 0) L 2 (1-2ε γ ,1) e C 2 t e M 2 0 4ε (α 2 -2α)t ≤ e M 0 α 2ε e -M 0 α(1-2ε γ ) 2ε θ ε,2 0 L 2 (1-2ε γ ,1) e C 2 t e M 2 0 4ε (α 2 -2α)t = θ ε,2 0 L 2 (1-2ε γ ,1) e C 2 t e -M 0 α 2ε (-2ε γ +(1-α 2 )M0t) , using that (recall that α > 0) e M 0 αx 2ε L ∞ (0,1) = e M 0 α 2ε and e -M 0 αx 2ε L ∞ (1-2ε γ ,1) = e -M 0 α(1-2ε γ ) 2ε
. The value α = ε 1-γ then leads to

θ ε,2 (•, t) L 2 (0,1) ≤ θ ε,2 0 L 2 (1-2ε γ ,1) e C 2 t e M0 e M 2 0 ε 1-2γ 4 t e -M 2 0 2ε γ t ≤ c 3 ε 1/2 e -M 2 0 2ε γ t (42)
for some constant c 3 > 0, using that γ ∈ (0, 1/2]. From inequalities (40) and (42) we deduce (36). The initial layer concentrated in the neighborhood of x = 1 vanishes exponentially fast with respect to t since the velocity transport is strictly positive.

3.1.2 Estimate of L ε P ε L 1 (0,t;L 2 (0,1))
Recall that L ε is the differential operator defined by (35). We have the following result proven in Appendice A.1.

Lemma 3.2. Assume hypotheses of Theorem 1.1. The function P ε defined by (32) satisfies

L ε (P ε ) = 5 i=1 I ε i in (0, 1) × (0, T ), with I ε 1 (x, t) = -εy 0 xx (x, t), I ε 2 (x, t) =ε 1/2 M xx (κ 1 , t) -M xx (a(t), t) w 2 2 W 0 w (w, t) + εM xx (κ 1 , t) w 2 2 W 1/2 w (w, t) -y 0 x (a(t) ± , t) , I ε 3 (x, t) = -y 0 t (1, t)e -M (1,t)z + M t (1, t) + M (1, t)M x (κ 2 , t) y 0 (1, t) -y 0 (a(t) ± , t) ze -M (1,t)z + M t (1, t) + M (1, t)M x (κ 2 , t) W 0 (τ M (1, t), t)ze -M (1,t)z -W 0 t (τ M (1, t), t) + τ M t (1, t)) - M (1, t) ε 1/2 W 0 w (τ M (1, t), t) e -M (1,t)z , I ε 4 (x, t) = -ε 1/2 W 0 wt (τ M (1, t), t) + τ M t (1, t) - M (1, t) ε 1/2 W 0 ww (τ M (1, t), t) ze -M (1,t)z -ε 1/2 M x (κ 2 , t)W 0 w (τ M (1, t), t)ze -M (1,t)z + ε 1/2 M t (1, t) + M (1, t)M x (κ 2 , t) W 0 w (τ M (1, t), t)z 2 e -M (1,t)z , I ε 5 (x, t) = ε 1/2 M t (1, t) + M (1, t)M x (κ 2 , t) C 1/2 (0, τ, t)ze -M (1,t)z -ε 1/2 C 1/2,t (0, τ, t) - C 1/2,τ (0, τ, t) ε 1/2 e -M (1,t)z , with inf {a(t), a(t) + ε 1/2 w} < κ 1 < sup{a(t), a(t) + ε 1/2 w}, 1 -εz < κ 2 < 1.
For the previous expansion, we get the following estimate for the L 1 (0, T ; L 2 (0, 1)) norm. Details of the proof are given in Appendix A.3. Lemma 3.3. Assume hypotheses of Theorem 1.1. Let P ε be the function defined by (32). There exists a constant c independent of ε such that

L ε P ε L 1 (0,t;L 2 (0,1)) ≤ c ε 1/2 , ∀t ∈ (0, T ]. (43)

Gronwall-Bellman type estimate

We now derive a priori estimates for the function z ε . Preliminary, since z ε is not vanishing at x = 0 and x = 1, we define

Z ε (x, t) := z ε (x, t) -z ε (0, t)f 0 ε (x) -z ε (1, t)f 1 ε (x), (x, t) ∈ Q T , (44)
where

f 0 ε ∈ C 2 ([0, 1], R + ) and satisfies f 0 ε (0) = 1 and f 0 ε (1) = 0, and f 1 ε (x) = f 0 ε (1 -x). The function Z ε solves the equation      L ε Z ε = L ε P ε -L ε (z ε (0, t)f 0 ε (x)) -L ε (z ε (1, t)f 1 ε (x)), (x, t) ∈ Q T , Z ε (0, t) = Z ε (1, t) = 0, t ∈ (0, T ), Z ε (x, 0) = -z ε (0, 0)f 0 ε (x) -z ε (1, 0)f 1 ε (x), x ∈ (0, 1), (45)
where z ε (0, 0) := lim t→0 + z ε (0, t), z ε (1, 0) := lim t→0 + z ε (1, t). Recall that L ε is defined by (35).

In order to derive an estimate for Z ε , we shall use the following Gronwall-Bellman type inequality (see [START_REF] Shao | Gronwall-Bellman type inequalities and their applications to fractional differential equations[END_REF]). Lemma 3.5. Let Z ε be the function defined by (44) and P ε the function defined by (32). There is a constant c independent of ε such that, for each t in [0, T ],

Lemma 3.4. Let I = [t 0 , T ] ⊂ R, k, b, p ∈ C(I, R + ). If ζ ∈ C(I, R + ) satisfies ζ(t) ≤ k(t) + t t0 b(s)ζ(s) ds + t t0 p(s)ζ γ (s) ds, t ∈ I, with 0 ≤ γ < 1, then for t ∈ I, ζ(t) ≤ A 1-γ (t) + (1 -γ) t t0 e (γ-1) t t 0 b(σ) dσ p(s) ds
Z ε (•, t) 2 L 2 (0,1) + 2ε t 0 Z ε x (•, s) 2 L 2 (0,1) ds ≤ |z ε (0, 0)| 2 f 0 ε 2 L 2 (0,1) + |z ε (1, 0)| 2 f 1 ε 2 L 2 (0,1) 1/2 +e Ct 2 t 0 1 0 L ε P ε (x, s) -L ε (f 0 ε (x)z ε (0, s) -L ε (z ε (1, s)f 1 ε (x))) 2 dx 1/2 ds 2 e Ct , (46)
with C := max

(x,t)∈Q T |M x (x, t)|.
Proof. Multiplying equation ( 45) by Z ε and integrating over (0, 1) × (0, t) gives

1 2 Z ε (•, t) 2 L 2 (0,1) + ε t 0 Z ε x (•, s) 2 L 2 (0,1) ds = 1 2 t 0 1 0 M x (x, s)Z ε (x, s) 2 dxds + 1 2 |z ε (0, 0)| 2 f 0 ε 2 L 2 (0,1) - 1 2 |z ε (1, 0)| 2 f 1 ε 2 L 2 (0,1) + t 0 1 0 L ε P ε (x, s) -L ε (z ε (0, s)f 0 ε (x)) -L ε (z ε (1, s)f 1 ε (x)) Z ε (x, s) dxds.
Applying the Cauchy-Schwarz inequality we obtain

Z ε (•, t) 2 L 2 (0,1) + 2ε t 0 Z ε x (•, s) 2 L 2 (0,1) ds ≤ C t 0 1 0 Z ε (x, s) 2 dxds + |z ε (0, 0)| 2 f 0 ε 2 L 2 (0,1) + |z ε (1, 0)| 2 f 1 ε 2 L 2 (0,1) + 2 t 0 1 0 L ε P ε (x, s) -L ε (z ε (0, s)f 0 ε (x)) -L ε (z ε (1, s)f 1 ε (x)) 2 dx 1/2 1 0 Z ε (x, s) 2 dx 1/2 ds, with C := max (x,s)∈Q T |M x (x, s)|. Applying Lemma 3.4 with t 0 = 0, b(s) = C, γ = 1 2 , A = k = |z ε (0, 0)| 2 f 0 ε 2 L 2 (0,1) + |z ε (1, 0)| 2 f 1 ε 2 L 2 (0,1) , ζ(t) = Z ε (•, t) 2 L 2 (0,1) + 2ε t 0 Z ε x (•, s) 2 L 2 (0,1) ds, p(s) = 2 1 0 L ε P ε (x, s) -L ε (z ε (0, s)f 0 ε (x) -L ε (z ε (1, s)f 1 ε (x)) 2 dx 1/2
, we obtain (46).

3.1.4 Estimate of z ε (x 0 , •) L 1 (0,s) , z ε t (x 0 , •) L 1 (0,s) and |z ε (x 0 , 0)|, x 0 = 0, 1
We now estimate each term of the right-hand side of (46).

Lemma 3.6. There is a constant c independent of ε such that, ∀s ∈ (0, T ],

z ε (0, •) L 1 (0,s) + z ε (1, •) L 1 (0,s) ≤ cε, (47)
z ε t (0, •) L 1 (0,s) + z ε t (1, •) L 1 (0,s) ≤ c, (48)
|z ε (x 0 , 0)| ≤ c, (49)
with z ε (x 0 , 0) := lim t→0 z ε (x 0 , t), x 0 = 0, 1.

Proof. Let s ∈ (0, T ] be arbitrary, and

w 0 (t) = -a(t) ε 1/2 , w 1 (t) = 1-a(t) ε 1/2 , t ∈ [0, T ].
We have from (33) that

z ε (0, t) = W 0 (w 0 (t), t) -v(0) -y 0 (X(0; 1, t)) + W 0 (τ M (1, t), t) -y 0 (0) e -M (1,t) ε + ε 1/2 W 1/2 (w 0 (t), t) -y 0 x (a(t) ± , t)w 0 (t) -ε 1/2 C 1/2 (0, τ, t) + 1 ε W 0 w (τ M (1, t), t) e -M (1,t) ε = W 0 (w 0 (t), t) -v(0) + ε 1/2 W 1/2 (w 0 (t), t) -y 0 x (a(t) ± , t)w 0 (t) + O(e -M (1,t) ε ).
We have

|W 0 (w 0 (t), t) -v(0)| = |y 0 (0) -v(0)| 2 erf c -w 0 (t) 2β(t) B(t) ≤ |y 0 (0) -v(0)|e - a 2 (t) 4εβ 2 (t)B(t) ,
then we deduce that W 0 (w 0 (t), t) -v(0) L 1 (0,s) ≤ cε.

By similar arguments to that used in the proof of Lemma 3.3 (point (e.2)), we show that W 1/2 (w 0 (t), t)y 0 x (a(t) ± , t)w 0 (t) L 1 (0,s) ≤ cε 1/2 allowing to conclude that z ε (0, •) L 1 (0,s) ≤ cε.

We also have

z ε (1, t) = W 0 (w 1 (t), t) -W 0 (τ M (1, t), t) + ε 1/2 W 1/2 (w 1 (t), t) -y 0 x (a(t) ± , t)w 1 (t) -ε 1/2 W 1/2 (τ M (1, t), t) -τ M (1, t)y 0 x (a(t) ± , t) . Writing W 0 (w 1 (t), t) -W 0 (τ M (1, t), t) = (W 0 (w 1 (t), t) -v(0)) -(W 0 (τ M (1, t), t) -v(0)
, and arguing as above we get that

W 0 (w 1 (t), t) -v(0) L 1 (0,s) + (W 0 (τ M (1, t), t) -v(0) L 1 (0,s) ≤ cε.
Similarly as above we have

W 1/2 (w 1 (t), t) -y 0 x (a(t) ± , t)w 1 (t) L 1 (0,s) ≤ cε 1/2 , W 1/2 (τ M (1, t)), t) -y 0 x (a(t) ± , t)τ M (1, t) L 1 (0,s) ≤ cε 1/2 , then conclude that z ε (1, •) L 1 (0,s) ≤ cε. Inequality (47) follows.
Concerning the second inequality, we have

z ε t (1, t) = W 0 t (w 1 (t), t) + w 1 (t)W 0 w (w 1 (t), t) -W 0 t (τ M (1, t), t) + τ M t (1, t) - M (1, t) ε 1/2 W 0 w (τ M (1, t), t) + ε 1/2 W 1/2 t (w 1 (t), t) + ε 1/2 w 1 (t) W 1/2 w (w 1 (t), t) -y 0 x (a(t) ± , t) + ε 1/2 M x (a(t), t)y 0 x (a(t) ± , t)w 1 (t) -ε 1/2 W 1/2 t (τ M (1, t), t) -ε 1/2 τ M t (1, t) - M (1, t) ε 1/2 W 1/2 w (τ M (1, t), t), t) -y 0 x (a(t) ± , t) -ε 1/2 M x (a(t), t) y 0 x a(t) ± , t τ M (1, t), t).
We have the following estimates, see the proof of Lemma 3.3 (points (c.3), (e.3), and (e.4), respectively),

W 0 t (τ M (1, t), t) + τ M t (1, t) - M (1, t) ε 1/2 W 0 w (τ M (1, t), t) L 1 (0,s) ≤ c, W 1/2 t (τ M (1, t), t) L 1 (0,s) ≤ cε 1/2 , W 1/2 w (τ M (1, t), t) -y 0 x (a(t) ± , t) L 1 (0,s) ≤ cε 1/2 .
Then clearly,

ε 1/2 τ M t (1, t) - M (1, t) ε 1/2 W 1/2 w (τ M (1, t), t), t) -y 0 x (a(t) ± , t) L 1 (0,s) ≤ cε 1/2 , ε 1/2 M x (a(t), t)y 0 x (a(t) ± , t)τ M (1, t), t) L 1 (0,s) + ε 1/2 M x (a(t), t)y 0 x (a(t) ± , t)w 1 (t) L 1 (0,s) ≤ c.
We have from (74)

W 0 t (w, t) + CW 0 w (w, t) = c --c + 2 1 2 √ π 2wα(t)γ(t) + w -2Cγ(t) γ(t) 3/2 e -w 2 4γ(t) , then taking w = w 1 (t) = 1-a(t) ε 1/2 and C = w 1 (t) = -a (t) ε 1/2 there holds that W 0 t (w, t) + CW 0 w (w, t) L 1 (0,t) = CO(ε 1/2 ) + O( ) = O(1), then W 0 t (w 1 (t), t) + w 1 (t)W 0 w (w 1 (t), t) L 1 (0,s) ≤ c.
We show, as above, that

ε 1/2 w 1 (t) W 1/2 w (w 1 (t), t) -y 0 x (a(t) ± , t) L 1 (0,s) ≤ cε 1/2 .
As in the proof of Lemma 3.3 (point (e.3)), we get ε 1/2 W 1/2 t (w 1 (t), t)

L 1 (0,s) ≤ cε. Thus z t (1, •) L 1 (0,s) ≤ c.
Similarly z t (0, •) L 1 (0,s) ≤ c, so that (48) follows. Eventually, we easily check that (49) holds. This ends the proof of the lemma.

End of the proof of Theorem 1.1

We are now in position to finish the proof of Theorem 1.1. It remains to choose the function f 0 ε ∈ C 2 ([0, 1]), satisfying f 0 ε (0) = 1 and f 0 ε (1) = 0 so as to minimize the terms in the right side of (46), asymptotically with respect to ε. We consider the functions

f 0 ε (x) = (1 -x)e -x ε , f 1 ε (x) = xe -1-x ε , x ∈ [0, 1], so that f i ε L 2 (0,1) ≤ c ε 1/2 and -εf i ε + M f i ε L 2 (0,1) ≤ c ε -1/2
, for i = 0, 1. It follows that, for x 0 = 0, 1, and for each t ∈ (0, T ],

f ε L 2 (0,1) z ε t (x 0 , •) L 1 (0,t) + |z ε (x 0 , x 0 )| ≤ cε 1/2 , -εf i ε + M f i ε L 2 (0,1) z ε (x 0 , •) L 1 (0,t) ≤ cε 1/2
. Coming back to Lemma 3.5, using the two previous estimates and Lemma 3.3, there holds that

Z ε (•, t) L 2 (0,1) + ε 1/2 Z ε x L 2 ((0,1)×(0,t)) ≤ cε 1/2 , ∀t ∈ [0, T ]. Now, since z ε (x, t) = Z ε (x, t) + f 0 ε (x)z ε (0, t) + f 1 ε (x)z ε (1, t), we deduce that z ε (•, t) L 2 (0,1) ≤ Z ε (•, t) L 2 (0,1) + |z ε (0, t)| f 0 ε L 2 (0,1) + |z ε (1, t)| f 1 ε L 2 (0,1) ≤ cε 1/2 + c |z ε (0, t)| + |z ε (1, t)| ε 1/2 , ∀t ∈ [0, T ]. ( 50
)
Writing z ε (x 0 , t) = z ε (x 0 , 0) + t 0 z ε t (x 0 , s) ds and using (49), there holds that

|z ε (x 0 , t)| ≤ |z ε (x 0 , 0)| + z ε t (x 0 , •) L 1 (0,t) ≤ c, ∀t ∈ [0, T ].
These estimate allow to deduce from (50) that z ε (•, t) L 2 (0,1) ≤ cε 1/2 for all t ∈ [0, T ]. Finally, since P ε -y ε = z ε -θ ε , using Lemma 3.1, we derive the estimate

P ε (•, t) -y ε (•, t) L 2 (0,1) ≤ c ε 1/2 + c e -M (1,0)ε γ ε + c ε 1/2 e -M 2 0 2ε γ t ≤ c ε 1/2 + c ε 1/2 e -M 2 0 2ε γ t , ∀t ∈ [0, T ].
This ends the proof of Theorem 1.1. 2

Some explicit examples

We explicit in this section some expansion of P ε for some particular case of the function M = M (x, t) then numerically evaluate the L ∞ (0, T ; L 2 (0, 1)) norm of the difference y ε -P ε .

4.1 The case M (x, t) = 1 + x, y 0 (x) ≡ y 0 in (0, 1) , v(t) ≡ v in (0, T )

The characteristic starting from (0, 0) is given by x -a(t) = 0 where the function a(t) solves

a (t) = M (a(t), t), t > 0, a(0) = 0,
and is given by a(t) = e t -1. T 0 such that a(T 0 ) = 1 is then given by T 0 = ln 2. Then,

α(t) := M x (a(t), t) = 1, β(t) := e t 0 α(s) ds = e t , B(t) := t 0 1 β(s) 2 ds = 1 -e -2t
2 , leading to 2β(t) B(t) = 2(e 2t -1). We obtain that y 0 defined by ( 7) is given by

y 0 (x, t) = y 0 , x + 1 > e t , v, x + 1 < e t ,
and

W 0 (w, t) = y 0 -v 2 erf w 2(e 2t -1) + y 0 + v 2 , W 1/2 (w, t) = 0, w = x + 1 -e t ε 1/2 .
p 0 defined by ( 24) is given by p 0 (x, t) = W 0 (w, t) while simply p 1/2 (x, t) = 0. Then, we get

C 0 (τ, t) = W 0 (2τ, t), C 1/2 (z, τ, t) = -zW 0 w (2τ, t) , z = 1 -x ε , τ = ln 2 -t ε 1/2 , implying P 0 (x, t) = W 0 (w, t) -W 0 (2τ, t)e -2z , P 1/2 (x, t) = -ze -2z W 0 w (2τ, t).
Consequently, our approximation is given by

P ε (x, t) = W 0 (w, t) -W 0 (2τ, t)e -2z -ε 1/2 ze -2z W 0 w (2τ, t), (x, t) ∈ Q T , w = x + 1 -e t ε 1/2 , z = 1 -x ε , τ = ln 2 -t ε 1/2 . ( 51
) Remark that ∂ w erf ( w g(t) ) = 2 √ π e -w 2 g(t) g(t) so that W 0 w (w, t) = y0-v √ π e - w 2 √ 2(e 2t -1) √ 2(e 2t -1)
. As an illustration, Figure 2 depicts the function P ε in Q T for ε = 10 -1 , 10 -2 , 5 × 10 -3 with T = 0.8 > T 0 . In particular, observe that P ε does not vanish at x = 0 and x = 1 but gets small as ε decreases. Table 1 provides some numerical values of the norm y ε -P ε L ∞ (0,T ;L 2 (0,1)) for several values of the parameter ε. The unknown solution y ε of ( 1) is obtained here using a numerical approximation. We obtain y ε -P ε L ∞ (0,T ;L 2 (0,1)) = O(ε 0.61 ) in full agreement with our estimate in Theorem 1.1. Table 1: y ε -P ε L ∞ (0,T ;L 2 (0,1)) w.r.t. ε for T = 0.8, y 0 = 1 and v = 0.

4.2

The case M (x, t) = 1 + t, y 0 (x) ≡ y 0 in (0, 1) , v(t) ≡ v in (0, T )

In this case, the function a(t) is given by a(t

) = t(t+2) 2 . T 0 such that a(T 0 ) = 1 is T 0 = √ 3 -1. Then, α(t) := M x (a(t), t) = 0, β(t) := e t 0 α(s) ds = 1 B(t) := t 0 1 β(s) 2 ds = t, leading to β(t) B(t) = √ t.
We obtain that y 0 defined by ( 7) is given by

y 0 (x, t) =      y 0 , x > t(t + 2) 2 , v, x < t(t + 2) 2 ,
and

W 0 (w, t) = y 0 -v 2 erf w 2 √ t + y 0 + v 2 , W 1/2 (w, t) = 0, w = 2x -t(t + 2) 2ε 1/2 .
p 0 defined by ( 24) is given by p 0 (x, t) = W 0 (w, t) while simply p 1/2 (x, t) = 0. Then, we get

C 0 (τ, t) = W 0 (2τ, t), C 1/2 (z, τ, t) = -zW 0 w (2τ, t) , z = 1 -x ε , τ = √ 3 -1 -t ε 1/2 , implying P 0 (x, t) = W 0 (w, t) -W 0 (1 + t)τ, t)e -(1+t)z , P 1/2 (x, t) = -ze -(1+t)z W 0 w ((1 + t)τ, t).
Consequently, our approximation is given

P ε (x, t) = W 0 (w, t) -W 0 ((1 + t)τ, t)e -(1+t)z -ε 1/2 ze -(1+t)z W 0 w ((1 + t)τ, t), w = 2x -t(t + 2) 2ε 1/2 , z = 1 -x ε , τ = √ 3 -1 -t ε 1/2 , ( 52
)
with here

W 0 w (w, t) = y0-v √ π e -w 2 2 √ t 2 √ t .

Concluding remarks

We have performed an asymptotic analysis of the solution of an advection-diffusion parametrized by a small coefficient in front of the diffusion term. The analysis takes into account the intersection of several singular layers and leads to an approximation in an ε 1/2 -neighborhood of the solution for the norm L ∞ (0, T ; L 2 (0, 1)). This analysis extends [START_REF]Asymptotic analysis of an advection-diffusion equation involving interacting boundary and internal layers[END_REF] to non constant, strictly positive transport velocity, leading to much technical arguments. In particular, the solution in the internal layer (following the characteristic starting form (0, 0)) solves a parabolic equation with non constant and unbounded coefficients related to the first derivatives of the velocity. The resulting solution is expressed in term of the heat kernel of the heat equation.

The ε 1/2 rate obtained for the final estimate is driven by our approximation W 0 in the internal layer, which depends on the choice non unique of the function g 0 . The choice [START_REF] Plaschko | Matched asymptotic approximations to solutions of a class of singular parabolic differential equations[END_REF] we made here is simple but notably involves the properties lim t→0 W 0 ( -a(t) ε 1/2 , t) = (y 0 (0) + v(0))/2 and lim x→0 W 0 ( x ε 1/2 , 0) = y 0 (0). This gap at the point (0, 0) generates an artificial boundary layer in the approximation of the solution along the line x = 0 above the characteristic x-a(t) = 0 in the neighborhood of t = 0. This boundary layer propagates inside the domain along the characteristic x -a(t) = 0 and affects the quality of the approximation. By analogy with [3, Section 3] devoted to the constant velocity case, one may consider the more involved choice

g ε 0 (w) :=    y 0 (0), w ≥ 0, v(0) + v(0) -y 0 (0) 2 e F w ε 1/2 , w < 0, ( 53
)
with

F := lim t→0 + a(t) β(t)B(t) = a (0) > 0 leading to W 0 ε = W 0 + U 0 ε with U 0 ε (w, t) = v(0) -y 0 (0) 2 e F w β(t)ε 1/2 + F 2 B(t) ε erf c w 2 B(t)β(t) + F B(t) ε 1/2 ,
and W 0 defined by [START_REF] Hubert | Vibration and coupling of continuous systems[END_REF]. The new function W 0 ε still satisfies the asymptotic conditions lim w→+∞ W 0 ε (w, t) = y 0 (0) and lim w→-∞ W 0 ε (w, t) = v(0) for all t ∈ (0, T ). Moreover, we now observe that

lim t→0 + W ε 0 -a(t) ε 1/2 , t = v(0), lim x→0 + W 0 ε x ε 1/2 , 0 = y 0 (0).
This allows to describe more accurately the discontinuity between the initial and Dirichlet conditions in the neighborhood of the point (0, 0). By analogy with [START_REF]Asymptotic analysis of an advection-diffusion equation involving interacting boundary and internal layers[END_REF], we may expect a O(ε 3/4 ) rate, but this remains to be checked. Similarly, as in [START_REF]Asymptotic analysis of an advection-diffusion equation involving interacting boundary and internal layers[END_REF], in order to get a better rate, one would needs to consider additional terms in the developments, notably the outer term y 1 defined by (2.2) not used here, internal layer terms W 1 , W 3/2 involving the heat kernel and boundary layer terms Y 1 and Y 3/2 solutions of ordinary differential equations with respect to the variable z. Whether or not we can perform the computations and get a priori estimates with these additional terms is an open issue.

In a futur work we plan to study the parabolic system describing the miscible displacement of compressible fluids in a porous medium. The displacement of one compressible fluid by another, completely miscible with the first, in a one-dimensional porous medium Ω = (0, 1), assuming for simplicity that the two fluids have the same compressibility factor z, is described by the following differential system, see for instance [START_REF] Bear | Introduction to modeling of transport phenomena in porous media[END_REF],

z∂ t p + ∂ x q = 0, q = - κ µ(c) ∂ x p, (x, t) ∈ Q T , (54)
φ∂ t c + q(x, t)∂ x c -∂ x (d(x, t)∂ x c)) = 0, (x, t) ∈ Q T , (55)
to which initial and boundary conditions are added. The unknowns of system (54), (55) are the functions p = p(x, t), c = c(x, t) and q = q(x, t); p is the pressure, c is the concentration of one of the two components of the fluid mixture, and q is the Darcy velocity. The function κ = κ(x) is the permeability of the medium, the constant φ is the porosity, µ = µ(c) is the concentration-dependent viscosity of the fluid mixture, d = d(x, t) is a small term representing the molecular diffusion and dispersion in the porous medium. The viscosity µ(c) is assumed to be determined by some mixing rule. See [START_REF] Amirat | Global weak solutions to equations of compressible miscible flows in porous media[END_REF] for the existence of a weak solution to the one-dimensional equations governing compressible flows of m miscible components in a porous medium. Our aim is to perform an asymptotic analysis of system (54), (55) when d is small and the concentration c satisfies Dirichlet boundary conditions on x = 0 and x = 1.

Acknowledgment. The authors thank the funding by the French government research program "Investissements d'Avenir" through the IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25).

A Appendix

A.1 Proof of Lemma 3.2

According to (30) and (31),

P ε (x, t) = P 0 (x, t) + ε 1/2 P 1/2 (x, t) =y 0 (x, t) + W 0 (w, t) -y 0 (a(t) ± , t) -y 0 (1, t) + W 0 (τ M (1, t), t) -y 0 (a(t) ± , t) e -M (1,t)z + ε 1/2 W 1/2 (w, t) -y 0 x (a(t) ± , t)w -ε 1/2 C 1/2 (0, τ, t) + W 0 w (τ M (1, t), t)z e -M (1,t)z .
Let us note here that the function P 0 + ε 1/2 P 1/2 belongs to C 1 ([0, 1] × (0, T ]), and that (P

0 -W 0 ) + ε 1/2 (P 1/2 -W 1/2 ) belongs to C 1 (Q T ).
Then, from the regularity assumptions on y 0 , v and M , (P

0 -W 0 ) + ε 1/2 (P 1/2 -W 1/2 ) belongs to H 2 (Q T ). We have L ε (W 0 ) = W 0 t (w, t) - a (t) ε 1/2 W 0 w (w, t) + M (x, t) ε 1/2 W 0 w (w, t) -W 0 ww (w, t) = W 0 t (w, t) + M (ε 1/2 w + a(t), t) -M (a(t), t) ε 1/2 W 0 w (w, t) -W 0 ww (w, t).
Using the Taylor expansion

M (ε 1/2 w + a(t), t) -M (a(t), t) = ε 1/2 wM x (a(t), t) + ε w 2 2 M xx (κ 1 , t), (56)
with inf {a(t), a(t) + ε 1/2 w} < κ 1 < sup{a(t), a(t) + ε 1/2 w}, and equation [START_REF] Deuring | L 2 -stability independent of diffusion for a finite elementfinite volume discretization of a linear convection-diffusion equation[END_REF] we get

L ε (W 0 ) = ε 1/2 w 2 2 M xx (κ 1 , t)W 0 w (w, t).
By similar calculations we get

L ε (W 1/2 ) = - w 2 2 M xx (a(t), t)W 0 w (w, t) + ε 1/2 w 2 2 M xx (κ 1 , t)W 1/2 w (w, t).
We observe that, when calculating L ε (P 0 + ε 1/2 P 1/2 ), it suffices to perform the calculation in Ω + ∪ Ω -, where

Ω + = {(x, t) ∈ Q T : x > a(t)} and Ω -= {(x, t) ∈ Q T : x < a(t)}. Straightforward calculations then give L ε y 0 (x, t) -y 0 (a(t) ± , t) = -εy 0 xx (x, t). (57)
Moreover,

L ε (y 0 x (a(t) ± , t)w) = w y 0 xt (a(t) ± , t) + a (t)y 0 xx (a(t) ± , t) - a (t) ε 1/2 y 0 x (a(t) ± , t) + M (x, t) ε 1/2 y 0 x (a(t) ± , t),
leading, using (56), to

L ε (y 0 x (a(t) ± , t)w) =w y 0 xt (a(t) ± , t) + M (a(t), t)y 0 xx (a(t) ± , t) + wM x (a(t), t) y 0 x (a(t) ± , t) + ε 1/2 w 2 2 M xx (κ 1 , t) y 0 x (a(t) ± , t).
Differentiating equation ( 6) 1 with respect to x yields

y 0 xt (x, t) + M x (x, t) y 0 x (x, t) + M (x, t) y 0 xx (x, t) = 0 for (x, t) ∈ Ω + ∪ Ω -.
Letting x → a(t) + and x → a(t) -, we get

y 0 xt (a(t) ± , t) + M x (a(t), t) y 0 x (a(t) ± , t) + M (a(t), t) y 0 xx (a(t) ± , t) = 0. Therefore L ε (y 0 x (a(t) ± , t)w) = ε 1/2 w 2 2 M xx (κ 1 , t) y 0 x (a(t) ± , t).
It results that

L ε W 0 (w, t) + ε 1/2 W 1/2 (w, t) -y 0 x (a(t) ± , t)w = ε w 2 2 M xx (κ 1 , t) W 1/2 w (w, t) -y 0 x (a(t) ± , t) + ε 1/2 w 2 2 W 0 w (w, t) M xx (κ 1 , t) -M xx (a(t), t) . (58)
We have

L ε y 0 (1, t)e -M (1,t)z = y 0 t (1, t)e -M (1,t)z -y 0 (1, t) M t (1, t) + M (1, t)M x (κ 2 , t) ze -M (1,t)z , L ε y 0 (a(t) ± , t)e -M (1,t)z = -y 0 (a(t) ± , t) M t (1, t) + M (1, t)M x (κ 2 , t) ze -M (1,t)z , L ε W 0 (τ M (1, t), t)e -M (1,t)z = W 0 t (τ M (1, t), t) + τ M t (1, t)) - M (1, t) ε 1/2 W 0 w (τ M (1, t), t) e -M (1,t)z -M t (1, t) + M (1, t)M x (κ 2 , t) W 0 (τ M (1, t), t)ze -M (1,t)z , with 1 -εz < κ 2 < 1. We deduce that L ε y 0 (1, t) + W 0 (τ M (1, t), t) -y 0 (a(t) ± , t) e -M (1,t)z = y 0 t (1, t)e -M (1,t)z -M t (1, t) + M (1, t)M x (κ 2 , t) y 0 (1, t) -y 0 (a(t) ± , t) ze -M (1,t)z + W 0 t (τ M (1, t), t) + τ M t (1, t) - M (1, t) ε 1/2 W 0 w (τ M (1, t), t) e -M (1,t)z -M t (1, t) + M (1, t)M x (κ 2 , t)W 0 (τ M (1, t), t) ze -M (1,t)z . (59)
We also have

L ε W 0 w (τ M (1, t), t)ze -M (1,t)z = W 0 wt (τ M (1, t), t) + τ M t (1, t)) - M (1, t ε 1/2 W 0 ww (τ M (1, t), t) ze -M (1,t)z + M x (κ 2 , t)W 0 w (τ M (1, t), t)ze -M (1,t)z -M t (1, t) + M (1, t)M x (κ 2 , t) W 0 w (τ M (1, t), t)z 2 e -M (1,t)z ,
and

L ε C 1/2 (0, τ, t)e -M (1,t)z = C 1/2,t (0, τ, t) - C 1/2,τ (0, τ, t) ε 1/2 e -M (1,t)z -C 1/2 (0, τ, t) M t (1, t) + M (1, t)M x (κ 2 , t) ze -M (1,t)z , (60)
where C 1/2,t (0, τ, t), C 1/2,τ (0, τ, t) are given by

C 1/2,t (0, τ, t) = W 1/2 t (τ M (1, t), t) + τ M t (1, t) W 1/2 w (τ M (1, t), t) -y 0 x (a(t) ± , t) + τ M (1, t)M x (a(t), t)y 0 x (a(t) ± , t), (61)
and

C 1/2,τ (0, τ, t) = M (1, t) W 1/2 w (τ M (1, t), t) -y 0 x (a(t) ± , t) . (62)
Collecting equalities (57)-( 60) we obtain the lemma.

A.2 An equality

Lemma A.1. Let K, β and B be the functions defined in Lemma 2.2. Let γ(τ ) := β 2 (τ )B(τ ) and

f (ξ, τ ) := 1 2 1 √ π 1 γ(τ ) 1/2 M xx (a(τ ), τ )ξ 2 e -ξ 2 4γ(τ ) , τ ∈ (0, t), ξ ∈ R. For all (w, t) ∈ R × (0, T ), t 0 R K(w, t; ξ, τ )f (ξ, τ ) dξdτ = 1 π e - w 2 4β 2 (t)B(t) β 2 (t)B 2 (t) B(t) t 0 M xx (a(τ ), τ ) B(τ ) B(τ )w 2 + 2β 2 (t)B 2 (t) -2β 2 (t)B(t)B(τ ) × × B(t) -B(τ )dτ. ( 63
)
Proof. We have, for all (w, t) ∈ R × (0, T ),

R K(w, t; ξ, τ )f (ξ, τ ) dξ = 1 √ π R e -v 2 f β(τ ) β(t) w + 2β(τ ) (B(t) -B(τ )) v, τ dv = 1 π 1 2γ(τ ) 1/2 M xx (a(τ ), τ ) R e -v 2 f β(τ ) β(t) w + 2β(τ ) B(t) -B(τ ) v, τ dv, with R e -v 2 f β(τ ) β(t) w + 2β(τ ) B(t) -B(τ ) v, τ dv = R e -v 2 e - β(τ ) β(t) w+2β(τ ) √ B(t)-B(τ ) v 2 4γ(τ ) β(τ ) β(t) w + 2β(τ ) B(t) -B(τ ) v 2 dv = 2e - w 2 4β(t) 2 B(t) β(τ ) B(τ ) B(τ )w 2 + 2β 2 (t)B 2 (t) -2β 2 (t)B(t)B(τ ) β 2 (τ )(B(t) -B(τ ))B(τ ) β 2 (t)B 2 (t) B(t) = 2e - w 2 4β(t) 2 B(t) β 2 (τ )B(τ ) B(τ )w 2 + 2β 2 (t)B 2 (t) -2β 2 (t)B(t)B(τ ) B(t) -B(τ ) β 2 (t)B 2 (t) B(t) ,
leading to (63).

A.3 Proof of Lemma 3.3

Let s ∈ (0, T ] be arbitrary. We estimate the L 1 (0, s; L 2 (0, 1)) norm of each term I ε i of the expression L ε (P ε ) = 5 i=1 I ε i , given by Lemma 3.2. In particular, we use that

z(x) n e -M (1,t)z(x) L ∞ (0,s;L 2 (Ω)) = O(ε 1/2 ) for n = 0, 1 and z(x) = 1 -x ε . (64)
Formula (63) will also be used several times. To avoid repetition, some points of the proof are shortened.

(a) Estimate of I ε 1 L 1 (0,s;L 2 (Ω)) . Clearly,

I ε 1 L 1 (0,s;L 2 (Ω)) ≤ cε. ( 65
) (b) Estimate of I ε 2 L 1 (0,s;L 2 (Ω)) . (b.1) Estimate of ε 1/2 M xx (κ 1 , t) -M xx (a(t), t) w 2 2 W 0 w (w, t) L 1 (0,s;L 2 (Ω))
. Let us denote

J 1 := ε s 0 1 0 1 2 x -a(t) ε 1/2 2 W 0 w x -a(t) ε 1/2 , t 2 dxdt.
Using the change of variable x-a(t) ε 1/2 = w, there holds that

J 1 = ε 3/2 s 0 1-a(t) ε 1/2 - a(t) ε 1/2 w 2 2 W 0 w (w, t) 2 dwdt ≤ ε 3/2 s 0 +∞ -∞ w 2 2 W 0 w (w, t) 2 dwdt,
provided that the last integral is finite. Using the explicit expression

W 0 w (w, t) = y 0 (0) -v(0) 2 2 √ π 1 2β(t) B(t) e - w 2
4β(t) 2 B(t) ,
and the change of variable

w 2β(t) √ B(t)
= v, there holds that

J 1 ≤ ε 3/2 s 0 +∞ -∞ 8β(t) 2 B(t) 3/2 c + -c - 2 1 √ π v 2 e -v 2 2 dvdt ≤ cε 3/2 .
We conclude that

ε 1/2 M xx (κ 1 , t) -M xx (a(t), t) w 2 2 W 0 w (w, t) L 1 (0,s;L 2 (Ω)) ≤ cε 3/4 . ( 66
) (b.2) Estimate of ε w 2 2 M xx (κ 1 , t) W 1/2 w (w, t) -y 0 x (a(t) ± , t) L 1 (0,s;L 2 (Ω))
. Let us denote

J 2 := ε 2 s 0 1 0 1 2 x -a(t) ε 1/2 2 W 1/2 w x -a(t) ε 1/2 , t -y 0 x (a(t) ± , t) 2 dxdt.
Using the change of variable x-a(t) ε 1/2 = w, there holds that

J 2 = ε 5/2 s 0 1-a(t) ε 1/2 - a(t) ε 1/2 w 2 2 W 1/2 w (w, t) -y 0 x (a(t) ± , t) 2 dwdt ≤ ε 5/2 s 0 +∞ -∞ w 4 W 1/2 w (w, t) -y 0 x (a(t) ± , t) 2 dwdt,
provided that the last integral is finite. We have from ( 23)

W 1/2 w (w, t) = g 1/2,w (w, t) + t 0 R K w (w, t; ξ, τ ) f 1/2 (ξ, τ ) -g 1/2 (ξ, τ ) dξdτ, w ∈ R, t ∈ (0, T ),
where g 1/2 is defined by ( 21) and g 1/2 is defined by [START_REF] Shao | Asymptotic analysis and domain decomposition for a singularly perturbed reaction-convectiondiffusion system with shock-interior layer interactions[END_REF]. Writing W 1/2 w (w, t) -y 0 x (a(t) ± , t) = W 1/2 w (w, t) -g 1/2,w (w, t) + g 1/2,w (w, t) -y 0 x (a(t) ± , t) , and noting that g 1/2,w (w, t) -y 0 x (a(t) ± , t) = 0 for |w| ≥ 2, we get

J 2 ≤ ε 5/2 s 0 +∞ -∞ w 4 W 1/2 w (w, t) -g 1/2,w (w, t) 2 dwdt + ε 5/2 .
Let us show that the integral

J 2 := s 0 +∞ -∞ w 4 W 1/2 w (w, t) -g 1/2,w (w, t)
2 dwdt is finite. We have

J 2 = s 0 +∞ -∞ w 4 t 0 +∞ -∞ K w (w, t; ξ, τ ) f 1/2 (ξ, τ ) -g 1/2 (ξ, τ ) dξdτ 2 dwdt.
Consider the integral (68)

J 1 2 := s 0 +∞ -∞ w 4 t 0 +∞ -∞ K w (w, t; ξ, τ )f 1/2 (ξ, τ ) dξdτ 2 dwdt. Differentiating iden- tity (63) in w gives t 0 +∞ -∞ K w (w, t; ξ, τ )f 1/2 (ξ, τ ) dξdτ = A 1 (w, t) + A 2 (w, t), with A 1 (w, t) := v(0) -y 0 (0) π (-w)e - w 2 4β 2 (t)B(t) 2β 4 (t)B 3 (t) B(t) × × t 0 M xx (a(τ ), τ ) B(τ ) B(τ )w 2 + 2β 2 (t)B 2 (t) -2β 2 (t)B(t)B(τ ) B(t) -B(τ )dτ, (67)
It results that

J 1 2 ≤ 2 s 0 +∞ -∞ w 4 A 1 (w, t) 2 + A 2 (w, t) 2 dwdt. (69)
Using the change of variable

w 2β(t) √ B(t) = v, there holds that s 0 +∞ -∞ w 4 A 1 (w, t) 2 dwdt = s 0 +∞ -∞ 32 v(0) -y 0 (0) π 2 v 6 e -2v 2 β(t)B 7/2 (t) × × t 0 M xx (a(τ ), τ ) B(τ ) B(τ )(4β 2 (t)B(t)v 2 ) + 2β 2 (t)B 2 (t) -2β 2 (t)B(t)B(τ ) × × B(t) -B(τ )dτ 2 dvdt,
then, using that B is strictly increasing, it results that

s 0 +∞ -∞ w 4 A 1 (w, t) 2 dwdt ≤ c s 0 +∞ -∞ β 3 (t)B 3/2 (t) t 0 |M xx (a(τ ), τ )| B(t) -B(τ )dτ 2 v 6 (v 2 + 1) 2 e -2v 2 dvdt, hence s 0 +∞ -∞ w 4 A 1 (w, t) 2 dwdt < ∞. Similarly, s 0 +∞ -∞ w 4 A 2 (w, t) 2 dwdt ≤ c s 0 +∞ -∞ β 3 (t)B
J 2 2 = s 0 +∞ -∞ w 4 t 0 +2 -2 K w (w, t; ξ, τ )g 1/2 (ξ, τ ) dξdτ 2 dwdt.
Explicitly,

K w (w, t; ξ, τ ) = -w β(t) -ξ β(τ ) 4 √ πβ(t)β(τ )(B(t) -B(τ )) 3/2 e - ( w β(t) - ξ β(τ ) ) 2 4(B(t)-B(τ )) , (w, ξ) ∈ R 2 , 0 ≤ τ < t ≤ T,
then, using the change of variable There are two positive constants c 1 and c 2 such that c 1 ≤ β(τ ) β(t) ≤ c 2 , 0 ≤ τ < t ≤ T. Let w 1 = 4 c1 , and w 2 = -4 c2 . We observe that:

ξ β(τ ) -w β(t) 2 √ (B(t)-B(τ )) = v gives +2 -2 K w (w, t; ξ, τ )g 1/2 (ξ, τ ) dξ = 1 4 √ πβ(t)β(τ )(B(t) -B(τ )) b(w,t,τ ) a(w,t,τ ) 2ve -v 2 g 1/2 β(τ ) β(t) w + 2β(τ ) B(t) -B(τ ) v, τ dv,
w ≥ w 1 =⇒ a(w, t, τ ) ≤ b(w, t, τ ) ≤ 0; w ≤ -w 2 =⇒ b(w, t, τ ) ≥ a(w, t, τ ) ≥ 0, then deduce that:

for w ≥ w 1 , b(w,t,τ ) a(w,t,τ ) 2|v|e -v 2 dv = e -b(w,t,τ ) 2 -e -a(w,t,τ ) 2 ; for w ≤ -w 2 , b(w,t,τ ) a(w,t,τ ) 2|v|e -v 2 dv = e -a(w,t,τ ) 2 -e -b(w,t,τ ) 2 .
For w ≥ w 1 , using the inequality -2 -β(τ ) β(t) w ≤ -β(τ ) β(t) w, we deduce that e -a(w,t,τ

) 2 ≤ e - (c 1 w) 2 4β 2 (t)(B(t)-B(τ ))
. Then there is a positive constant c 3 such that

1 β(t)β(τ )(B(t) -B(τ )) e -a(w,t,τ ) 2 ≤ c 3 e - (c 1 w) 2 8β 2 (t)(B(t)-B(τ )) . (70)
We also have, for

w ≥ w 1 , 2 -β(τ ) β(t) w 2 ≥ 1 2 β(τ ) β(t) w 2 ≥ 1 2 c 1 w 2 , then we deduce that 1 β(t)β(τ )(B(t) -B(τ )) e -b(w,t,τ ) 2 ≤ c 3 e - (c 1 w) 2 8β 2 (t)(B(t)-B(τ )) . (71)
Inequalities ( 70) and (71) hold similarly for w ∈ (-∞, w 2 ) so that, for w ∈ (-∞, w 2 ) ∪ (w 1 , ∞),

t 0 +2 -2 K w (w, t; ξ, τ )g 1/2 (ξ, τ ) dξdτ ≤ c t 0 e - (c 1 w) 2 8β 2 (t)(B(t)-B(τ )) dτ.
Using the Cauchy-Schwarz inequality, there holds that

t 0 +2 -2 K w (w, t; ξ, τ )g 1/2 (ξ, τ ) dξdτ 2 ≤ cT t 0 e - (c 1 w) 2 4β 2 (t)(B(t)-B(τ )) dτ, then +∞ -∞ |w| 4 t 0 +2 -2 K w (w, t; ξ, τ )g 1/2 (ξ, τ ) dξdτ 2 dw ≤ cT t 0 +∞ -∞ |w| 4 e - (c 1 w) 2 4β 2 (t)(B(t)-B(τ )) dτ dw.
Using the change of variable

w 2β(t) √ (B(t)-B(τ )) = v, we get that t 0 +∞ -∞ |w| 4 e - (c 1 w) 2 4β 2 (t)(B(t)-B(τ )) dτ dw = t 0 2 5 β(t) 5 (B(t) -B(τ )) 5/2 +∞ -∞ |v| 4 e -c 2 1 v 2 dτ dv, then clearly, J 2 2 < ∞.
Using the Young inequality, we deduce that J 2 < ∞, then conclude that J 2 ≤ cε 5/2 , and

ε w 2 2 M xx (κ 1 , t) W 1/2 w (w, t) -y 0 x (a(t) ± , t) L 1 (0,s;L 2 (Ω)) ≤ cε 5/4 . (72)
It results from ( 66) and (72) that

I ε 2 L 1 (0,s;L 2 (Ω)) ≤ cε 3/4 . (73)
(c) Estimate of I ε 3 L 1 (0,s;L 2 (Ω)) . (c.1) Using (64) there holds that

y 0 t (1, t)e -M (1,t)z L 1 (0,s;L 2 (Ω)) ≤ cε 1/2 , M (1, t)M x (κ 2 , t) y 0 (1, t) -y 0 (a(t) ± , t) ze -M (1,t)z L 1 (0,s;L 2 (Ω)) ≤ cε 1/2 .
(c.2) From the explicit expression of W 0 , see [START_REF] Hubert | Vibration and coupling of continuous systems[END_REF], we get the uniform bound |W 0 (τ M (1, t), t)| ≤ max{|v(0)|, |y 0 (0)|} then, using (64), there holds that

M t (1, t) + M (1, t)M x (κ 2 , t) W 0 (τ M (1, t), t)ze -M (1,t)z L 2 (0,s;L 2 (Ω)) ≤ cε 1/2 . (c.3) Let γ(t) := β 2 (t)B(t)
. By a direct calculation we have, for all C ∈ R,

W 0 t (w, t) + CW 0 w (w, t) = v(0) -y 0 (0) 2 1 2 √ π 2wα(t)γ(t) + w -2Cγ(t) γ(t) 3/2 e -w 2 4γ(t) . ( 74
) Taking w = τ M (1, t) = T1-t ε 1/2 M (1, t) and C = τ M t (1, t) -M (1,t) ε 1/2 = T1-t ε 1/2 M t (1, t) -M (1,t) ε 1/2
, we get

W 0 t (w, t) + CW 0 w (w, t) L 1 (0,s) = CO(ε 1/2 ) + O( ) = O(1), then W 0 t (τ M (1, t), t) + τ M t (1, t)) - M (1, t) ε 1/2 W 0 w (τ M (1, t), t) e -M (1,t)z L 1 (0,s;L 2 (Ω)) ≤ cε 1/2 .

29

From the previous estimate and the estimates in (c.1) and (c.2) we deduce that

I ε 3 L 1 (0,s;L 2 (Ω)) ≤ cε 1/2 . ( 75
) (d) Estimate of I ε 4 L 1 (0,s;L 2 (Ω)) . (d.1) A direct calculation gives W 0 wt (w, t) + CW 0 ww (w, t) = v(0) -y 0 (0) 2 1 4 √ π -2w 2 α(t)γ(t) -w 2 + 2wCγ(t) + 4γ(t) 2 α(t) + 2γ(t) γ(t) 5/2 e -w 2 4γ(t) . (76)
Taking

w = τ M (1, t) = T1-t ε 1/2 M (1, t) and C = τ M t (1, t) -M (1,t) ε 1/2 = T1-t ε 1/2 M t (1, t) -M (1,t) ε 1/2
, we get

W 0 wt (τ M (1, t), t) + CW 0 ww (τ M (1, t), t) L 1 (0,s) = O(ε 1/2 ) + CO( ) = O(ε 1/2 ), hence ε 1/2 W 0 wt (τ M (1, t), t) + τ M t (1, t) - M (1, t) ε 1/2 W 0 ww (τ M (1, t), t) ze -M (1,t)z L 1 (0,s;L 2 (Ω)) ≤ cε 3/2 . (d.2) From W 0 w (w, t) = y0(0)-v(0) 2 1 √ π 1 γ(t) 1/2 e -w 2 4γ(t) with w = τ M (1, t) = T1-t ε 1/2 M (1, t), we get W 0 w (τ M (1, t), t) L 1 (0,s) = O(ε 1/2 ), ( 77
) then ε 1/2 M x (κ 2 , t)W 0 w (τ M (1, t), t)ze -M (1,t)z L 1 (0,s;L 2 (Ω)) ≤ cε 3/2 , ε 1/2 M t (1, t) + M (1, t)M x (κ 2 , t) W 0 w (τ M (1, t), t)z 2 e -M (1,t)z L 1 (0,s;L 2 (Ω)) ≤ cε 3/2 .
From the two last estimates and that in (d.1) we deduce that

I ε 4 L 1 (0,s;L 2 (Ω)) ≤ cε 3/2 . ( 78
)
(e) Estimate of I ε 5 L 1 (0,s;L 2 (Ω)) . (e.1) Using (64) we deduce that τ M (1, t)M x (a(t), t)y 0

x (a(t) ± , t) L 1 (0,s) = O(ε -1/2 ).

(e.2) Estimate of C 1/2 (0, τ, t) L 1 (0,s) . We use that 63) there holds that

W 1/2 (w, t) -y 0 x (a(t) ± , t)w = W 1/2 (w, t) -g 1/2 (w, t) + g 1/2 (w, t) -y 0 x (a(t) ± , t)w , with g 1/2 (w, t) -y 0 x (a(t) ± , t)w = 0 for |w| ≥ 2. Taking w = T1-t ε 1/2 M (1, t) in (
J 1 3 (t) := t 0 +∞ -∞ K T 1 -t ε 1/2 M (1, t), t; ξ, σ f 1/2 (ξ, σ) dξdσ = v(0) -y 0 (0) π e - T 1 -t ε 1/2 M (1,t) 2 4β 2 (t)B(t) β 2 (t)B 2 (t) B(t) × × t 0 M xx (a(σ), σ) B(σ) B(σ) T 1 -t ε 1/2 M (1, t) 2 + 2β 2 (t)B 2 (t) -2β 2 (t)B(t)B(σ) × × B(t) -B(σ)dσ.
Consequently,

J 1 3 (t) ≤ cT e - T 1 -t ε 1/2 M (1,t) 2 4β(t) 2 B(t) B(t) T 1 -t ε 1/2 M (1, t) 2 + 1 ≤ cT e - T 1 -t ε 1/2 M 1 2 2β 2 2 B 2
, where we used the bounds 0 < M 1 ≤ M (x, t), B(t) ≤ B 2 and β(t) ≤ β 2 for all (x, t) ∈ Q T . Then

T1 0 J 1 3 (t) 2 dt ≤ cT 2 T1 0 e - T 1 -t ε 1/2 M 1 2 β 2 2 B 2 dt.
Using the change of variable T1-t ε 1/2 = τ , there holds that

T1 0 J 1 3 (t) 2 dt ≤ cT 2 ε 1/2 +∞ 0 e - (τ M 1 ) 2 β 2 2 B 2 dτ ≤ cT 2 ε 1/2 .
When s ≥ T 1 , we derive similarly the estimate

s T1 J 1 3 (t) 2 dt ≤ cT 2 ε 1/2 , then get s 0 J 1 3 (t) 2 dt ≤ cT 2 ε 1/2 . Let J 2 3 (t) := t 0 +∞ -∞ K T1-t ε 1/2 M (1, t), t; ξ, σ g 1/2 (ξ, σ) dξdσ. We show similarly that s 0 J 2 3 (t) 2 dt ≤ cT 2 ε 1/2 ,
then derive the estimate

W 1/2 (τ M (1, t), t) -τ M (1, t)y 0 x (a(t) ± , t) L 1 (0,s) ≤ cT 2 ε 1/2 , then C 1/2 (0, τ, t) L 1 (0,s) = O(ε 1/2 ),
and the first term of I 5 satisfies the estimate

ε 1/2 M t (1, t) + M (1, t)M x (κ 2 , t) C 1/2 (0, τ, t)ze -M (1,t)z L 1 (0,s;L 2 (Ω)) ≤ cε 3/2 . (e.3) Estimate of' W 1/2 t (τ M (1, t), t) L 1 (0,s) . Let us denote F 1 (w, t) := v(0)-y0(0) π e - w 2 4β 2 (t)B(t) β 2 (t)B 2 (t) √ B(t) , F 2 (w, t) := t 0 M xx (a(σ), σ) B(σ) B(σ)w 2 + 2β 2 (t)B 2 (t) -2β 2 (t)B(t)B(σ) B(t) -B(σ)dσ,
and F (w, t) := F 1 (w, t)F 2 (w, t). We have from (63), for all (w, t) ∈ R × (0, T ),

t 0 R K(w, t; ξ, σ)f 1/2 (ξ, σ) dξdσ = F (w, t).
Differentiating the previous equality in t yields

t 0 R K t (w, t; ξ, σ)f 1/2 (ξ, σ) dξdσ + f 1/2 (w, t) = F t (w, t).
We have

F 1,t (w, t) = v(0) -y 0 (0) π e - w 2 4β 2 (t)B(t) w 2 2β (t)B(t) + β(t)B (t) 4β 5 (t)B 9/2 (t) - 2β (t)B(t) + 5/2B (t)β(t) β 3 (t)B 7/2 (t) , F 2,t (w, t) = t 0 M xx (a(σ), σ) B(σ) × × 4β(t)β (t)B 2 (t) + 4B(t)B (t)β 2 (t) -4β(t)β (t)B(t)B(σ) -2β 2 (t)B (t)B(σ) B(t) -B(σ) dσ + t 0 M xx (a(σ), σ) B(σ) B(σ)w 2 + 2β 2 (t)B 2 (t) -2β 2 (t)B(t)B(σ) 1 2 (B(t) -B(σ)) -1/2 B (t)dσ. Differentiating (23) in t gives W 1/2 t (w, t) = g 1/2,t (w, t) + t 0 R K t (w, t; ξ, σ) f 1/2 (ξ, σ) -g 1/2 (ξ, σ) dξdσ + f 1/2 (w, t) -g 1/2 (w, t). Then W 1/2 t (w, t) = g 1/2,t (w, t) + F t (w, t) - t 0 R K t (w, t; ξ, σ)g 1/2 (ξ, σ) dξdσ -g 1/2 (w, t). (79)
Arguing as in (e.2), using the bounds 0 < M = v, there holds that 2 B 2 where we used the bounds 0 < M 1 ≤ M (x, t), B(t) ≤ B 2 and β(t) ≤ β 2 for all (x, t) ∈ Q T . Then ≤ c.

(e.5) Estimate of ε 1/2 C 1/2,t (0, τ, t)e -M (1,t)z L 1 (0,s;L 2 (Ω)) . Using (64) and the estimates in (e.3) and (e.4), we deduce that ε 1/2 C 1/2,t (0, τ, t)e -M (1,t)z L 1 (0,s;L 2 (Ω))

≤ cε 1/2 .

(e.6) We have from (e.4) that W 1/2 w (τ M (1, t), t) -y 0 x (a(t) ± , t)

L 1 (0,s)

≤ cε 1/2 , then, using (64) we deduce that C 1/2,τ (0, τ, t)e -M (1,t)z L 1 (0,s;L 2 (Ω))

≤ cε 1/2 .

From the previous estimate, the last estimate in (e.2), and the estimate in (e.5) we deduce that

I ε 5 L 1 (0,s;L 2 (Ω)) ≤ cε 1/2 . ( 85
)
Collecting estimates (65), (73), (75), (78), and (85) we get the desired result.

Figure 1 :

 1 Figure 1: Internal (red) and boundary (blue) layer zones for y ε for a constant velocity M (x, t) = M > 0 (left) and non constant velocity M (x, t) > 0 (right).

Figure 2 :

 2 Figure 2: Function P ε given by (51) in Q T for ε = 10 -1 , 10 -2 , 5 × 10 -3 .

A 2 (β 2 0 M

 220 w, t) := v(0) -y 0 (t)B 2 (t) B(t) t xx (a(τ ), τ ) B(τ )B(τ ) B(t) -B(τ )dτ.

. 2 K

 2 Since g 1/2 (ξ, τ ) is uniformly bounded, we have+2 -w (w, t; ξ, τ )g 1/2 (ξ, τ ) dξ ≤ c β(t)β(τ )(B(t) -B(τ )) b(w,t,τ ) a(w,t,τ ) 2|v|e -v 2dv.

2 4 ( 1 β

 41 B(t)-B(σ)) , with R (w, t; ξ, σ) = -(σ) 4π(B(t) -B(σ)) (t) -B(σ)) -(t) -B(σ)) 2

2 B 2

 2 dτ ≤ cε 1/2 , using the change of variable

  3/2 (t)

						t	|M xx (a(τ ), τ )| B(t) -B(τ )dτ	2	v 6 e -2v 2	dvdt,
						0
	hence	s 0	+∞ -∞ w 4 A 2 (w, t)	2 dwdt < ∞. It results from (69) that the integral J 1 2 is finite.
	Consider the integral J 2 2 :=	s 0	+∞ -∞ w 4 t 0	+∞ -∞ K w (w, t; ξ, τ )g 1/2 (ξ, τ ) dξdτ	2	dwdt. Since g 1/2 (ξ, τ ) = 0
	for |w| ≥ 2 we have

  1 ≤ M (x, t), B(t) ≤ B 2 and β(t) ≤ β 2 for all (x, t) ∈ Q T , we find that|F 1 (w, t)| ≤ ce (w 2 + 1), |F 2 (w, t)| ≤ c(w 2 + 1), |F 2,t (w, t)| ≤ c + c(w 2 + 1)Using inequality (80) and the change of variable T1-t ε 1/2 = τ , there holds that When s ≥ T 1 , we derive similarly the estimate M (1, t), t dt ≤ cε 1/2 , then get K t (w, t; ξ, σ)g 1/2 (ξ, σ) dξdσ. We have K t (w, t; ξ, σ) = R (w, t; ξ, σ) e -

				-w 2 4β 2 2 B 2 , |F 1,t (w, t)| ≤ ce -w 2 4β 2 2 B 2 -w 2 4β 2 2 B 2 (w 2 + 1) ≤ ce	-w 2 4β 2 2 B 2 ,	(80)
		0	T1	F t	T 1 -t ε 1/2 M (1, t), t dt ≤ cε 1/2	0	+∞	e	-	(τ M 1 ) 2 β 2 2 B 2 dτ ≤ cε 1/2 .
					0	s T1 F t ε 1/2 s T1-t F t T 1 -t ε 1/2 M (1, t), t dt ≤ cε 1/2 .	(82)
	Consider the term	t 0 R ( w β(t)	-	ξ β(σ) )

t 0 (B(t) -B(σ)) -1/2 B (t) dσ = c + c(w 2 + 1) t 0 (B(t) -B(σ)) -1/2 B (σ) B (t) B (σ) dσ ≤ c(w 2 + 1).

It results that

|F t (w, t)| ≤ ce

for some constant c > c. Taking w = T1-t ε 1/2 M (1, t) in (79) we have, for ε small enough,

W 1/2 t T 1 -t ε 1/2 M (1, t), t = F t T 1 -t ε 1/2 M (1, t), t -t 0 R K t T 1 -t ε 1/2 M (

1, t), t; ξ, σ g 1/2 (ξ, τ ) dξdσ. (81)

  K t (w, t; ξ, σ)g 1/2 (ξ, σ) dξ -v 2 A 1 (t, σ)(2v 2 + 1) + A 2 (t, σ)|vw| dv.Using the inequalities e -v 2 (2v2 + 1) ≤ ce -v 2 2 , e -v 2 |v| ≤ ce -v 2 2, we haveK t (w, t; ξ, σ)g 1/2 (ξ, σ) dξ ≤ c A 1 (t, σ) + A 2 (t, σ)|w| (t, σ) + A 2 (t, σ)|w| A 3 (t, σ)e -a(w,t,σ) 2 M (1, t), t; ξ, σ g 1/2 (ξ, σ) dξ ≤ c T 1 -t ε 1/2 M (1, t) + 1 e -When s ≥ T 1 , we derive similarly the estimate M (1, t), t; ξ, σ g 1/2 (ξ, σ) dξdσ ≤ cε 1/2 , then getK t (w, t; ξ, σ)g 1/2 (ξ, σ) dξdσ dt ≤ cε 1/2 . ± , t) = W 1/2 w (w, t) -g 1/2,w (w, t) + g 1/2,w (w, t) -y 0 x (a(t) ± , t) , then taking w = T1-t ε 1/2 M t (1, t) we have, for ε small enough, M t (1, t), t; ξ, τ f 1/2 (ξ, τ ) dξdτ = A 1 T 1 -t ε 1/2 M t (1, t), t + A 2 T 1 -t ε 1/2 M t (1, t), t ,where A 1 (w, t) and A 2 (w, t) are defined by (67) and (68). We have |A 1 (w, t)| + |A 2 (w, t)| ≤ ce

	+2 -2 = 2β(σ) B(t) -B(σ) × g 1/2 w, t; β(σ) β(t) w + 2β(σ) B(t) -B(σ) v, σ dv, b(w,t,σ) a(w,t,σ) e -v 2 R w, t; β(σ) β(t) w + 2β(σ) B(t) -B(σ) v, σ × with a(w, t, σ) = -2-β(σ) β(t) w 2β(σ) √ B(t)-B(σ) and b(w, t, σ) = 2-β(σ) β(t) w 2β(σ) √ B(t)-B(σ) . We have 2β(σ) B(t) -B(σ) R w, t; β(σ) β(t) w + 2β(σ) B(t) -B(σ) v, σ = 1 √ π   B (t)(2v 2 -1) 2(B(t) -B(σ)) -vwβ (t) β 2 (t) B(t) -B(σ)   , then +2 -2 K t (w, t; ξ, σ)g 1/2 (ξ, σ) dξ = 1 √ π b(w,t,σ) a(w,t,σ) e -v 2 g 1/2 β(σ) β(t) w + 2β(σ) (B(t) -B(σ)) v, σ × b(w,t,σ) a(w,t,σ) e -v 2   B (t)(2v 2 + 1) 2(B(t) -B(σ)) + |vw|β (t) β 2 (t) B(t) -B(σ)   dv. Let us denote A 1 (t, σ) := B (t) 2(B(t) -B(σ)) , A 2 (t, σ) := β (t) β 2 (t) so that +2 -2 K t (w, t; ξ, σ)g 1/2 (ξ, σ) dξ ≤ c b(w,t,σ) a(w,t,σ) +2 -2 b(w,t,σ) a(w,t,σ) e -v 2 2 dv Integrating in t and using the change of variable T1-t ε 1/2 = τ , there holds that T1 0 t 0 +2 -2 K t T 1 -t ε 1/2 M (1, t), t; ξ, σ g 1/2 (ξ, σ) dξdσdt ≤ cε 1/2 T 1 ε 1/2 0 (τ M (1, t) + 1) t 0 e -a(τ M (1,t),t,σ) 2 4 dσdτ ≤ cε 1/2 +∞ 0 (τ M (1, t) + 1) t 0 e -a(τ M (1,t),t,σ) 2 4 dσdτ ≤ cε 1/2 . s T1 +2 -2 K t T1-t ε 1/2 s 0 t 0 R We conclude that W 1/2 t 1/2 We have from (23) W 1/2 Writing W 1/2 w (w, t) -y 0 x (a(t) W 1/2 w T 1 -t ε 1/2 M t (1, t), t -y 0 x (a(t) ± , t) = W 1/2 w T 1 -t ε 1/2 M t (1, t), t -g 1/2,w T 1 -t ε 1/2 M t (1, t), t . Then We have from (b.2) ≤ c A 1 Then +2 -2 K t T 1 -t ε 1/2 a T 1 -t ε 1/2 M (1,t),t,σ 2 4 . Integrating in σ, we have t 0 +2 -2 K t T 1 -t ε 1/2 M (1, t), t; ξ, σ g 1/2 (ξ, σ) dξdσ ≤ c T 1 -t 0 ε 1/2 M (1, t) + 1 t e -a T 1 -t ε 1/2 M (1,t),t,σ 2 4 dσ. (83) J 1 4 (t) : = t 0 +∞ -∞ T 1 -t K w ε 1/2 -w 2 2β 2
	W 1/2
	(84)

  B (t)(2v 2 -1) 2(B(t) -B(σ)) -vwβ (t) β 2 (t) B(t) -B(σ)   dv.

Using that g 1/2 (ξ, τ ) is uniformly bounded, there holds that

+2 -2 K t (w, t; ξ, σ)g 1/2 (ξ, σ) dξ ≤ c B(t) -B(σ) , A 3 (t, σ) := b(w, t, σ) -a(w, t, σ), e 4 .

Let us take w = T1-t ε 1/2 M (1, t) in the previous inequality. For ε small enough and t ≤ T 1 we have a(w, t, σ) < 0, then

A 1 (t, σ) + A 2 (t, σ)|w| A 3 (t, σ)e -a(w,t,σ) 2 2 ≤ c |w| + 1 e -a(w,t,σ)

2 4

.

(τ M (1, t), t)

L 1 (0,s) ≤ cε 1/2 . (e.4) Estimate of' τ M t (1, t) W w (τ M (1, t), t) -y 0 x (a(t) ± , t) L 1 (0,s)

. Let us denote

J 4 (t) := W 1/2 w T 1 -t ε 1/2 M t (1, t), t -y 0 x (a(t) ± , t). w (w, t) = g 1/2,w (w, t) + t 0 R K w (w, t; ξ, τ ) f 1/2 (ξ, τ ) -g 1/2 (ξ, τ ) dξdτ, w ∈ R, t ∈ (0, T ). w T 1 -t ε 1/2 M t (1, t), t -y 0 x (a(t) ± , t) = t 0 R K w T 1 -t ε 1/2 M

t (1, t), t; ξ, τ f 1/2 (ξ, τ ) -g 1/2 (ξ, τ ) dξdτ.

  T1-t ε 1/2 = τ . When s ≥ T 1 , we derive similarly the estimate M (1, t), t; ξ, σ g 1/2 (ξ, σ) dξdσ. We show similarly that s 0 J 2 4 (t) dt ≤ cε 1/2 , then derive the estimatesW 1/2 w (τ M (1, t), t) -y 0x (a(t) ± , t)

				s T1 J 1 4 (t)	2 dt ≤ cε 1/2 ,
	then get	s 0 J 1 4 (t) dt ≤ cε 1/2 . Let J 2 4 (t) =	t 0	+∞ -∞ K T1-t ε 1/2

L 1 (0,s) ≤ cε 1/2 , τ M t (1, t) W 1/2 w (τ M (1, t), t) -y 0 x (a(t) ± , t)

L 1 (0,s)


domain: spi phys

text:
Cyprien Brulon
email: cyprien.brulon@onera.fr

Baptiste Fix

Clément Verlhac

Patrick Bouchon

Printed terahertz metasurfaces for multispectral imaging by thermo-conversion

In this work, we investigate the fabrication and the application of ultra-thin terahertz metasurfaces as thermal converters for indirect terahertz imaging. The microstructures fabricated by an ultrasonically driven printing process (Microplotter) are used to improve the THz to IR conversion efficiency and tune the spectral or polarisation selectivity. The resulting conversion membranes show optical and thermal responses which are consistent with numerical simulations establishing reliable rules to design such membrane at various wavelengths. This work paves the way for a low-cost solution of multispectral terahertz imaging with a standard infrared camera.

I -INTRODUCTION

Terahertz (THz) imaging is an active field of research and finds numerous applications such as non-destructive industrial testing, security, chemical identification, medical imaging [START_REF] Dhillon | The 2017 terahertz science and technology roadmap[END_REF]. However the high cost, the low sensitivity and spectral selectivity of terahertz detectors constitute a serious obstacle for a wide range of applications. In this perspective, terahertz imaging by thermal conversion promises to be a low-cost alternative technique to conventional terahertz detectors but still suffers from low sensitivity (P min (1THz)=50µW/cm 2 ) in comparison with micro bolometers and CMOS based detectors [START_REF] Lewis | A review of terahertz detectors[END_REF].

In this work, we investigate the use of ultra-thin terahertz metasurfaces as a thermal converter and a way to improve the THz to IR conversion efficiency and tune the spectral/polarisation selectivity (see Fig 1 .). We demonstrate the fabrication of periodical arrays of metallic resonators such as split ring (SRR) on an ultra-thin substrate of PET (3µm) by an ultrasonically driven printing process (Microplotter). This modular method has been used for the rapid and low-cost fabrication of resonant structures in the 50GHz-1THz range [START_REF] Salmon | Rapid prototyping of flexible terahertz metasurfaces using a microplotter[END_REF]. The resulting conversion membranes show optical and thermal responses which are consistent with numerical simulations establishing reliable rules to design such membrane at various wavelengths, paving the way for a low-cost solution of multispectral terahertz imaging with a standard infrared camera.

Fig 1: Concept of terahertz imaging by thermoconversion.

The thz wave is absorbed by the membrane which then heats up. This heating is then measurable by a thermal camera. The advantage of such a detection system is that it transfers all the technological complexity to the infrared camera.

II -FABRICATION METHOD

As a low-cost manufacturing process an inkjet printing technique based on ultrasonic excitation was used instead of optical lithography to fabricate the membrane (Microplotter II) [START_REF] Larson | Controlled deposition of picoliter amounts of fluid using an ultrasonically driven micropipette[END_REF]. Unlike conventional printers, the metallic patterns are drawn in a pencil like fashion through the use of a glass tip (see Fig 2a .). After being filled with ink containing silver nanoparticles (JS-A291 -Metalon) by capillary action, the stretched glass capillary is brought close to a substrate. Thanks to a piezoelectric element, we measure the variation of the resonance frequency of the capillary once it is in contact. The annealing of the ink by thermal conduction (100 to 250 degrees for 5min) sets the quality of the metal (see Fig 2b .) with a compromise to be found: the more the ink is heated the better the conductivity will be which increases the quality factor of the terahertz resonance at the cost of a lower infrared emissivity of the metal. The source being naturally polarized, we have two ways of illuminating our sample either in transverse electric or in transverse magnetic knowing that only one of these polarizations gives rise to the appearance of a resonance for SRR type structures [START_REF] Pendry | Magnetism from conductors and enhanced nonlinear phenomena[END_REF] (see Fig 3b .). The electromagnetic response is simulated using COMSOL. To observe the thermal response of the membrane we must then replace the terahertz sensor with an infrared camera (FLIR SC700) showing that the heating is only present when there is resonance of the membrane (see Fig 3c .)

IV -DISCUSSION

The resulting conversion membrane is then naturally spectrally and polarisation sensitive with selectivity depending on the quality factor of the plasmonic resonators. For SRR-type structures, this quality factor (Q=λ/Δλ) is close to 20, which gives a selectivity of only a few gigahertz at 75 GHz. Other types of resonators such as guided mode resonators [START_REF] Verlhac | Towards perfect metallic behavior in optical resonant nanostructures[END_REF] or coupled nano fabry perot [START_REF] Fix | High-quality-factor double Fabry-Perot plasmonic nanoresonator[END_REF] can improve this selectivity of the conversion membrane by increasing the quality factor of the resonance. Moreover, this concept is transposable to any terahertz frequency by adapting the geometrical parameters of the resonators.

Finally, the membrane can be used in a terahertz imaging scheme with a near state-ofthe-art detectability threshold of 50µW/cm 2 and a maximum acquisition frequency at high terahertz power (>50mW/cm 2 ) of 1Hz. For the proof of concept, a metallic triangle hidden behind a 5 mm thick cardboard is scanned and imaged with a surprisingly good resolution of two to three times the wavelength used (2mm), which is very promising for low-cost nondestructive testing applications (see Fig 4 .).

Fig 2 :

 2 Fig 2: a. Diagram of the printing method. b. Scanning electron microscopy image of the ink after annealing at 350° for 5min.

Fig 3 :

 3 Fig 3: a. Unit cell of the printed terahertz metasurface consisting of split ring resonators sized to resonate at 75 GHz in TM bias (D=3000µm, L=2000µm, w=80µm and g=50µm). b. Experimental spectra measured by the sub-terahertz specrtoscopy bench in TE and TM polarization. c. Corresponding infrared image in TE and TM polarization.

Fig 4 :

 4 Fig 4: Therahertz image of a metal object hidden behind a 5 mm thick cardboard. This object is resolved at 2λ thanks to the ability of the terahertz waves to pass through the cardboard.

V -CONCLUSION

To conclude, we have demonstrated the use of inkjet printed terahertz metasurfaces as converters for terahertz imaging by thermoconversion. The fabricated microstructures have provided spectral and polarisation sensitivity in a terahertz imaging scheme. This work thus paves the way for a low-cost multispectral terahertz imaging solution based on a standard infrared camera.