Keywords: Problèmes de diffraction, milieux non bornées, modèles asymptotiques, impédances géneralisées, forte conductivité

Nous analysons la précision et le cractère bien posé des problèmes aux limites d'impédances généralisées dans le contexte des problèmes de diffraction par des milieux non bornés fortement absorbants. Nous nous restreignons dans ce premier travail au problème scalaire (polarization E). Comparé aux travaux précédents, le caractère non borné du milieu introduit une difficltée supplémentaire puisque les arguments classiques de compacité ne peuvent être utilisés. Notre analyse repose sur l'utilisation d'identités de Rellich et des estimations a priori de type L 2 sur la solution du problème. Nous discutons également l'approximation numérique des conditions obtenues et testons leur ordre de précision.

Introduction

Time harmonic wave scattering from rough layers is an important problem in science and engineering, as it describes for instance scattering of electromagnetic waves from the ground when one models the earth as a rough stratified medium. In such a model, the moisture of soil causes absorption of the electromagnetic wave inside the ground, and thus naturally leads to a scattering problem for a rough absorbing layer. Since waves inside the absorbing part of the medium decay exponentially with respect to the distance to the layer's boundary, a lot of research has been carried out how to replace the wave scattering problem inside the absorbing layer by some easily handable absorbing boundary condition on the interface in between the absorbing layer and free space [START_REF] Artola | Diffraction d'une onde électromagnétique par un obstacle borné à permittivité et perméabilité élevées[END_REF][START_REF] Duruflé | Higher order generalized impedance boundary conditions in electromagnetic scattering problems[END_REF][START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF][START_REF] Rytov | Calcul du skin-effet par la méthode des perturbations[END_REF][START_REF] Senior | Approximate boundary conditions in electromagnetics[END_REF]. The aim of such a boundary condition is to set up an approximate scattering problem merely in the complement of the absorbing object, while still guaranteeing a reliable error bound on the solution of the approximate problem. This error bound depends on what we call the order of the boundary condition as well as on the magnitude of the absorption inside the layer. Indeed, we treat the magnitude of absorption as a parameter and expand the acoustic field in a power series with respect to the inverse of this parameter. Approximate boundary conditions are built after truncation of this series, the order of obtained conditions then corresponds to the truncation index. Truncation at order zero simply leads to a Dirichlet boundary condition, which is naturally the formal limit condition as the absorption tends to infinity; truncation at order one leads to a (usual) impedance boundary condition. This is the reason why we call the conditions arising from truncation at higher order generalized impedance boundary conditions (GIBC).

In this paper, we analyse GIBC for rough absorbing layers up to order three and shall restrict ourselves to the scalar problem (which corresponds in 2-D to the E-polarization of electromagnetic waves). While the construction of such conditions is rather analogous to the case of a bounded absorbing inhomogeneity, the error analysis is more complicated. This is already obvious when one considers merely existence and uniqueness of solution for the approximate problems. For instance, for the case of a bounded obstacle, existence of solution for the time harmonic exterior Dirichlet or impedance scattering problem is known for a long time [START_REF] Nédélec | Acoustic and electromagnetic equations[END_REF]. For the rough surface scattering problem with a Dirichlet boundary condition, corresponding results have only been achieved during the last decade, firstly by using integral equation approach [START_REF] Chandler-Wilde | Acoustic scattering by mildly rough surfaces in three dimensions[END_REF][START_REF] Chandler-Wilde | Scattering by rough surfaces: the Dirichlet problem for the Helmholtz equation in a non-locally perturbed half-plane[END_REF][START_REF] Zhang | Integral equation methods for scattering by infinite rough surfaces[END_REF], and more recently, by using a variational approach in [START_REF] Chandler-Wilde | Existence, uniqueness, and variational methods for scattering by unbounded rough surfaces[END_REF][START_REF] Chandler-Wilde | The mathematics of scattering by unbounded, rough, inhomogeneous layers[END_REF]. For scattering from rough infinite layers we refer to recent results in [START_REF] Lechleiter | Scattering of acoustic waves from rough layers[END_REF]. For the variational theory on rough surface scattering, Rellich identities have been shown to be particularly useful since they provide a-priori bounds on a solution to the scattering problem, thereby establishing existence and uniqueness of solution via an inf-sup condition. Such a-priori bounds are also important for proving existence of solution to rough layer scattering problems involving higher order impedance boundary conditions. Moreover, they permit to construct a bounded L 2 solution operator. By definition, the solution operator maps Dirichlet boundary data on the rough surface to the radiating solution of the Helmholtz equation taking this boundary data. We show that this operator has a bounded extension from square integrable functions on the interface into the space of square integrable functions in a layer of finite height above the interface.

The important role of this L 2 solution operator in our analysis is to replace compactness arguments present in earlier rigorous error analysis of generalized impedance boundary conditions. Since compact embeddings of Sobolev spaces do not hold in our unbounded setting, we cannot use such types of compactness arguments which are always present in earlier works on generalized impedance boundary conditions for bounded objects [START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF]. The L 2 solution operator does part of this job. The other part is mainly done by a Rellich identity for radiating solutions of the Helmholtz equation over a rough layer. Through our Rellich identity we are able to prove existence and uniqueness of solution to the rough layer scattering problem subject to generalized impedance boundary conditions up to order three. Further, we show optimal error bounds for solutions to the approximate scattering problems involving our generalized impedance boundary conditions compared to the solution of the original scattering problem in the absorbing layer.

The structure of this paper is as follows. The first section is dedicated to the presentation of the mathematical setting of the scattering problem from unbounded rough surfaces and the introduction of used notation. Section 3 serves as a brief review of the main steps in deriving generalized impedance boundary conditions and needed extensions to the case of rough surfaces. Afterwards, we analyse lower order Neumann-to-Dirichlet conditions in Section 4. This cases enables us to present the analysis in a simpler framework. The case of more complicated higher order condition is analysed in Section 4. We indicate in particular how "stablized" conditions can be treated in a similar way as standard impedance boundary conditions. Finally, in Appendix A we provide some auxiliary existence and regularity results on rough surface scattering, some of which might have an interest in its own.

Settings of the problem and notation

Let us start with a brief description of the geometrical setting and our notation, such that we can afterwards present the problem mathematical setting in detail. Points in the Euclidean space R m (m = 2, 3) are denoted by x = (x 1 , x 2 , . . . , x m ) ⊤ and sometimes it is convenient to write x = (x, x m ) ⊤ , that is, x are the first m -1 coordinates of x ∈ R m . In this work, m = 2 or 3, although all arguments also work in higher dimension. By R m ± := {x ∈ R m , x m ≷ 0} we denote the upper and lower half space of R m and the plane in between

R m ± is called Γ 0 = {x ∈ R m , x m = 0}. More generally, Γ a = {x = (x 1 , x 2 , . . . , x m ) ⊤ ∈ R m , x m = a}, a ∈ R. The half space above and below Γ a is denoted by U ± a := {x ∈ R m , x m ≷ a}. The domain Ω := {x ∈ R m , -a < x m < a} is partitioned into two parts Ω ± := {x ∈ Ω, x n ≷ f (x)} by the interface Γ := {x ∈ R m , f (x) = x m },
given by a twice continuously differentiable function f : R m-1 → R, -a < f < a. For simplicity, we also introduce Ω

R := {x ∈ Ω, x 2 1 + . . . x 2 m-1 < R 2 }.
The exterior unit normal field to Ω and Ω R is called ν; on Γ we choose the unit normal field ν to point downwards into Ω -. The boundary

of Ω R is C R := ∂Ω R . We split C R = C + R ∪ C - R with C - R := {x ∈ Γ -a , x 2 1 + . . . x 2 m-1 < R 2
} and define M R := C R ∩ Ω. We refer to Figure 1 for a sketch of the geometry of the problem. where u| ± Γ is the limit taken from Ω ± . The L 2 inner product on Γ is denoted as •, • L 2 (Γ) . Let us also remark that the square root of the complex unit -i in the fourth quadrant is denoted by α = (1i)/ √ 2. All fields in this paper are time harmonic, that is, their time dependence is exp(-iωt) for frequency ω > 0, and this time dependence will always be suppressed. The wave number k is defined as k = ω/c with c > 0 the speed of sound in vacuum. Now we turn to the mathematical formulation of the rough layer problem. We describe the medium of propagation by a refractive index function n 2 , which is assumed to be a real valued function in Ω + and of the special form

x m ν ν ν Γ Γ a Γ -a Ω + Ω -
n 2 = 1 + i k 2 ε 2
in Ω -for ε > 0.

The (small) parameter ε will control the magnitude of absorption in Ω -throughout this paper.

The refractive index hence jumps across the interface; concerning smoothness, we always require that n 2 ∈ C 0,1 (Ω + ) and for some statements we even require more smoothness of n 2 Ω+ . Also, we always suppose that Re(n 2 ) ≥ c 0 > 0 and Im(n 2 ) ≥ 0 in Ω + . In the upper half space U + a we suppose n 2 to equal a real positive constant

n 2 + with Re(n 2 + ) ≥ c 0 such that n 2 Γa = n 2 + , while, obviously, n 2 equals n 2 -:= 1 + i/(k 2 ε 2 ) in U - -a .
Important for existence of solutions to rough layer scattering problems is the further assumption that ∂n 2 /∂x m ≤ 0 in Ω + . All these assumptions are supposed to hold throughout the paper. We remark that we could also deal with a refractive index whose real part varies in Ω -, but for simplicity do not consider this case here.

The total field due to a local time harmonic source g supported in Ω + satisfies the Helmholtz

equation ∆u + k 2 n 2 u = g in R m , (1)
subject to the additional assumption that u and its normal derivative be continuous over the interface Γ where the index of refraction n 2 jumps, and a radiation condition in U ± ±a . We note that the solution u to the above problem will depend on ε and denote u = u ε .

Let us now introduce the radiation condition imposed on u ε . As shown in [START_REF] Chandler-Wilde | Existence, uniqueness, and variational methods for scattering by unbounded rough surfaces[END_REF], the Fourier transform

F : L 2 (R m-1 ) → L 2 (R m-1 ), F φ(ξ) = (2π) -(m-1)/2 R m-1 e -ix•ξ φ(x) dx, ξ ∈ R m-1 , for φ ∈ L 2 (R m-1 ) ∩ L 1 (R m-
u ε (x) = (2π) -(m-1)/2 R m-1 exp i (x m -a) k 2 n 2 + -ξ 2 + x • ξ F u| Γa (ξ) dξ, (2)
for x ∈ U + a . Note that this representation is a superposition of upwards propagating plane waves. Computing the normal derivative of the latter expression reveals that the Dirichlet-to-Neumann operator

T + n 2 + : H 1/2 (Γ a ) → H -1/2 (Γ a ) is given by (T + n 2 + φ)(x) = i(2π) -(m-1)/2 R m-1 k 2 n 2 + -ξ 2 exp (ix • ξ) F φ(ξ) dξ and it is shown in [3] that T + n 2 + is bounded from H 1/2 (Γ a ) to H -1/
u ε,3 T + n 2 + (u ε,3 ) ds = Im i R m-1 k 2 -ξ 2 F ( u ε,3 Γa ) 2 ds = |ξ|<k k 2 -ξ 2 F ( u ε,3 Γa ) 2 ds ≥ 0 (3)
and a similar computation shows that -Re Γa u ε,3 T + n 2 + (u ε,3 ) ds ≥ 0; again, the same inequalities

holds for T - n 2 - .
The two Dirichlet-to-Neumann operators allow to frame the acoustic scattering problem under investigation variationally in the domain Ω. Due to the homogeneous jump conditions [u] Γ = [∂u/∂ν] Γ = 0 it makes sense to seek a solution u ε ∈ H 1 (Ω). Twice (formally) applying Green's identity in Ω ± shows that

Ω ∇u ε • ∇v -k 2 n 2 u ε v dx - Γa v T + n 2 + (u ε ) ds - Γ-a v T - n 2 - (u ε ) ds = - Ω gv dx (4)
for all v ∈ H 1 (Ω).

We now introduce the Neumann-to-Dirichlet operator D ε on Γ which maps φ ∈ H -1/2 (Γ) to the Dirichlet boundary values v ε on Γ of the unique solution to

∆v ε + k 2 + i ε 2 v ε = g in Ω -, ∂v ε /∂ν = T - n 2 - (v ε ) on Γ -a , ∂v ε /∂ ν = -φ on Γ.
Due to absorption, the map D ε : φ → v ε ∈ H 1/2 (Γ) is well defined and bounded. Of course, the restriction u ε | Ω+ solves

∆u ε + k 2 u ε = g in Ω + , ∂u ε /∂ν = T + n 2 + (u ε ) on Γ a , u ε + D ε (∂u ε /∂ν) = 0 on Γ. ( 5
)
Roughly speaking, the idea of a generalized impedance boundary condition is now to construct a (formal) expansion of D ε in terms of ε and to obtain explicitly computable boundary operators replacing D ε in [START_REF] Chandler-Wilde | The mathematics of scattering by unbounded, rough, inhomogeneous layers[END_REF]. These approximations to D ε will introduce a certain error which we will show to be bounded in terms of powers of ε, the actual power depending on the truncation of the asymptotic development of D ε . To conclude this brief outlook, as an approximation to the restriction u ε | Ω+ we are going to study problems of the following form,

∆u ε,p + k 2 n 2 u ε,p = g in Ω + , ∂u ε,p /∂ν = T + n 2 + (u ε,p ) on Γ,
together with u ε,p + D ε,p (∂u ε,p /∂ν) = 0 on Γ for certain Neumann-to-Dirichlet operators D ε,p of order p ∈ {0, 1, 2, 3}, which will be constructed and analysed in detail.

INRIA 3 Formal construction of GIBC

As a brief introduction to the construction of generalized impedance boundary conditions, we recall some basic ideas, definitions and computations from [START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF]. The formal computation of these conditions for rough layers is the same as for bounded absorbing inclusion and we can skip most of the computations. Let Ω δ -:= {x ∈ Ω -, dist(x, Γ) < δ} for δ > 0 small enough such that each point x ∈ Ω δ -has a unique representation x = x Γ + qν, with x Γ ∈ Γ and q > 0. Since we assumed that the function f which defines Γ is of class C 2 , such a choice δ is always possible. Recall that the unit normal ν on Γ was defined to point into Ω -. The parameter δ is fixed from now on and we also fix a cut off function

χ ∈ C ∞ (Ω -) such that χ = 1 in Ω δ/2 -, χ = 0 in Ω δ -as well as 0 ≤ χ ≤ 1, |∇χ| ≤ C and |∆χ| ≤ C in Ω -.
The starting point for the construction of generalized impedance boundary conditions is the assumption that the exact solution u ε of the scattering problem (4) with parameter ε > 0 can be written as

u ε (x) = u 0 + (x) + εu 1 + (x) + ε 2 u 2 + (x) + . . . , x ∈ Ω + , (6)
for functions u ℓ + : Ω + → C and

χu ε (x) = u 0 -(x Γ , q/ε) + εu 1 -(x Γ , q/ε) + ε 2 u 2 -(x Γ , q/ε) + . . . , x = x Γ + qν ∈ Ω δ -, (7)
with functions u ℓ -:

Γ×R + → C such that lim η→∞ u ℓ -(x Γ , η) = 0 for x Γ ∈ Γ.
From this expansion we will in the sequel construct boundary conditions which allow to truncate problem (4) at Γ, introducing an error that can be controlled in powers of ε. Following [START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF], we set

ũε -(x Γ , η) = u 0 -(x Γ , η) + εu 1 -(x Γ , η) + ε 2 u 2 -(x Γ , η) + . . . , x Γ ∈ Γ, η > 0.
Since u ε solves the Helmholtz equation, all functions u ℓ + also need to satisfy this equation. Moreover, starting from the Helmholtz equation ( 1) one can compute a differential equation for ũε -of the form

- ∂ 2 ∂η 2 -i ũε -= 8 ℓ=1 ε ℓ A ℓ ũε - on Γ × R + ,
where the A ℓ are differential operators in (x Γ , η) independent of ε; see [START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF] for details. Substituting the latter equation into [START_REF] Duruflé | Higher order generalized impedance boundary conditions in electromagnetic scattering problems[END_REF] shows that

- ∂ 2 ∂η 2 -i u p -= 8 ℓ=1 A ℓ u p-ℓ - in Γ × R + , (8)
where u p-ℓ -= 0 for p -ℓ < 0. Let us think of this equation as a family of second order ordinary differential equations in η with parameter x Γ . Coupling of (8) with the expansion (6) of u ε in Ω + yields a boundary condition at η = 0 and together with the decay condition for u ℓ -we can solve [START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF]. In more detail, equating terms in ( 6) and ( 7) which share the same powers of ε either offers the possibility of coupling the Dirichlet traces,

u p + (x Γ ) = u p + (x Γ , 0), for x Γ ∈ Γ, (9)
or coupling of the normal derivatives,

∂u p-1 + ∂ν (x Γ ) = ∂u p + ∂η (x Γ , 0), for x Γ ∈ Γ. ( 10
)
We concentrate in this paper on the second option, which leads to Neumann-to-Dirichlet impedance boundary conditions, and which is somehow more natural due to the shift of the index. The first option ( 9) results in Dirichlet-to-Neumann impedance boundary conditions which will not be considered here; the interested reader is referred to [START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF] where formulas for Dirichlet-to-Neumann impedance boundary conditions are provided which can be transferred to the rough surface context by the techniques presented in this paper. The differential equations in η in (8) together with [START_REF] Lechleiter | Scattering of acoustic waves from rough layers[END_REF] and the decay condition u p -(x Γ , η) → 0 as η → ∞ can be explicitly solved in the form

u p -(x Γ , η) = P p
xΓ (η)e -αη , η > 0, p ∈ N 0 , with a polynomial P p xΓ of degree p depending on ∂u ℓ -/∂ν, ℓ ∈ {0, 1, . . . , p -1}. There holds a recurrence relation of order 8 for the P p xΓ , which we will however not give here, but only give the explicit form of the u p -for p ∈ {0, 1, 2, 3}. Before citing the result of this computation in [8, Eq. (4.26)-(4.29)], we remark that we denote Gauss and mean curvature on Γ by G and H, respectively; we also assume familiarity of the reader with standard surface differential operators such as the surface gradient ∇ Γ and the surface Laplacian ∆ Γ , see, e.g, [START_REF] Nédélec | Acoustic and electromagnetic equations[END_REF] for details. Using the abbreviation α = (1i)/ √ 2 there holds

u 0 -(x Γ , η) =0, (11)
u 1 -(x Γ , η) = - 1 α ∂u 0 + ∂ν (x Γ )e -αη , (12)
u 2 -(x Γ , η) = - 1 α ∂u 1 + ∂ν (x Γ ) + H α 2 ∂u 0 + ∂ν (x Γ ) + η H α ∂u 0 + ∂ν (x Γ ) e -αη , (13)
u 3 -(x Γ , η) = - 1 α ∂u 2 + ∂ν (x Γ ) + H α 2 ∂u 1 + ∂ν (x Γ ) - 1 2α 3 (3H 2 -G + k 2 ) ∂u 0 + ∂ν (x Γ ) (14)
- 1 2α 3 ∆ Γ ∂u 0 + ∂ν + η H α ∂u 1 + ∂ν (x Γ ) - 1 2α 2 (∆ Γ + 3H 2 -G + k 2 ) ∂u 0 + ∂ν (x Γ ) +η 2 1 2α (G -3H 2 ) ∂u 0 + ∂ν (x Γ ) e -αη .
In the following lemma, we show that the functions u p + for p ∈ {0, 1, 2, 3} are well defined in H 1 (Ω + ) if the right hand side g of the Helmholtz equation ( 1) is smooth enough. Lemma 1. Set g 0 = g ∈ L 2 (Ω + ) and g ℓ = 0 for ℓ = 1, 2, . . . and let p ∈ N = {1, 2, . . . }. For g ∈ H p-1 (Ω + ), n 2 ∈ C p-1,1 (Ω + ) and Γ of class C r,1 , the functions u ℓ + are well defined in H 1 (Ω + ) for ℓ = 0, 1, . . . , p through the recursion

∆u ℓ + + k 2 n 2 u ℓ + = g ℓ in Ω + , ∂u ℓ + ∂ν = T + n 2 + (u ℓ + ) on Γ a , u ℓ + = u ℓ -(•, 0) on Γ, ( 15
)
where u ℓ -(•, 0) is determined from - ∂ 2 ∂η 2 -i u ℓ -= 8 ℓ=1 A ℓ u ℓ-l - in Γ × R + , lim η→∞ u ℓ -(x Γ , η) = 0 on Γ,
together with the coupling condition [START_REF] Lechleiter | Scattering of acoustic waves from rough layers[END_REF]. Additionally, the bounds

u 0 + H 2 (Ω+) ≤ C g L 2 (Ω+) , ∂u 0 + /∂ν H 1/2 (Γ) ≤ C g L 2 (Ω+)
as well as

u ℓ + H p-ℓ (Ω+) ≤ C g H p-1 (Ω+) and ∂u ℓ + /∂ν H p-(1+2ℓ)/2 (Γ) ≤ C g H p-1 (Ω+) hold as long as p -(1 + 2ℓ)/2 > 0. INRIA Proof.
For ℓ = 0 we know from ( 11) that u 0 -(•, 0) = 0. Hence, u 0 + solves a homogeneous Dirichlet problem for the Helmholtz equation with right hand side g 0 = g. From Theorem 11 we conclude that u 0 + ∈ H 2 (Ω) for g ∈ L 2 (Ω + ); thus the normal derivative ∂u 0 + /∂ν belongs to H 1/2 (Γ) and the bounds ∂u 0 + /∂ν H 1/2 (Ω+) ≤ C u 0 + H 2 (Ω+) ≤ C g L 2 (Ω+) hold. More generally, for g ∈ H p-1 (Ω + ) and n 2 as well as Γ smooth enough we have that u 0 + ∈ H p+1 (Ω); thus the normal derivative ∂u 0 + /∂ν belongs to H p-1/2 (Γ) and the bounds

∂u 0 + /∂ν H p-1/2 (Γ) ≤ C u 0 + H p+1 (Ω+) ≤ C(p) g H p-1 (Ω+)
hold for p = 1, 2, . . . . In consequence, the Dirichlet trace u 1 -(•, 0) given in ( 12) is well defined in H p-1/2 (Γ) for g in H p (Ω + ) and we can solve the Dirichlet problem [START_REF] Zhang | Integral equation methods for scattering by infinite rough surfaces[END_REF] in H p (Ω + ) according to Theorem 11, with

∂u 1 + /∂ν H p-3/2 (Γ) ≤ C ∂u 1 + /∂ν H p (Γ) ≤ C g H p-1 (Ω+)
, as long as p -3/2 > 0. Existence of u ℓ + for ℓ = 2, . . . , p follows from an induction argument. Note that existence and boundedness of

u ℓ + in H 1 (Ω + ) and ∂u ℓ + /∂ν ∈ H 1/2 (Γ) requires g ∈ H p-1 (Ω + ) with p -(1 + 2ℓ)/2 > 0, as well as n 2 ∈ C p-1,1 (Ω + ) and Γ of class C p,1 .
Now we can derive a first version of generalized impedance boundary conditions. Except for the case p = 3 these will be the ones for which we prove convergence of optimal order in Section 4. For p = 3, sophisticated further manipulations of the impedance boundary condition are necessary to derive convergence of optimal order. Therefore we devote the entire section 5 to the analysis of the case p = 3.

The basic idea behind an impedance boundary condition of order p is truncation of the expansion [START_REF] Chandler-Wilde | Scattering by rough surfaces: the Dirichlet problem for the Helmholtz equation in a non-locally perturbed half-plane[END_REF] of

u ε at ℓ = p, ũε,p = p ℓ=0 ε ℓ u ℓ + in Ω + . (16)
In view of ( 15) we get that ũε,p = p ℓ=0 11)-( 14) in terms of the normal derivatives ∂u ℓ + /∂ν, we can plug in these formulas into the last equation to obtain for ℓ = 0 that ũε,0 = 0 and for ℓ = 1 that ũε,1 + ε

ε ℓ u ℓ -(•, 0) on Γ. Since u ℓ -(•, 0) is explicitly given in (
1 α ∂u 0 + ∂ν (x Γ ) = 0, that is, ũε,1 + ε 1 α ∂ ũε,1 ∂ν (x Γ ) = ε 2 1 α ∂ ũε,1 ∂ν .
Generally speaking, we find a tangential differential operator D ε,p acting on Γ such that

ũε,p + D ε,p ∂ ũε,p ∂ν = ε p+1 r ε,p , p ∈ {0, 1, 2, 3}, (17)
for functions r ε,p which are explicitly given in [8, Eq. (4.34)] for p ∈ {0, 1, 2, 3} as

r ε,0 = 0, r ε,1 = 1 α ∂u 1 + ∂ν , r ε,2 = 1 α ∂u 2 + ∂ν -iH ∂ ∂ν u 1 + + εu 2 + , r ε,3 = 1 α ∂u 3 + ∂ν -iH ∂ ∂ν u 1 + + εu 2 + - 1 2 ∆ Γ + 3H 2 -G + k 2 ∂ ∂ν (u 1 + + εu 2 + + ε 2 u 3 + ) . ( 18
)
The differential operators D ε,p are given by [8, Eq. (3.5)-(3.8)]

D ε,0 = 0, D ε,1 = ε α , D ε,2 = ε α -iε 2 H, D ε,3 = ε α -iε 2 H -ε 3 α 2 (∆ Γ + 3H 2 -G + k 2 ).

RR n°6841

Neglecting the small right hand side in (17), we define an approximation u ε,p of u ε , solution of (4), by

∆u ε,p + k 2 n 2 u ε,p = g in Ω + , ∂u ε,p ∂ν = T + n 2 + (u ε,p ) on Γ a , u ε,p + D ε,p u ε,p = 0 on Γ, (19)
for p ∈ {0, 1, 2}. For p = 3, u ε,3 will be defined in a more refined way in Section 5. We note that this u ε,p is defined through a boundary value problem merely posed in Ω + . Well-posedness of this problem, as well as the approximation of u ε by u ε,p is subject of the following sections. However, at least morally we already note from (17) that the condition of order p will introduce an error u ε,p -u ε which is O(ε p ).

Impedance Boundary Conditions of Lower Order

Our aim in this section is to provide basic tools for the convergence analysis of generalized impedance boundary conditions leading to proofs of convergence of u ε,p to u ε with optimal order p + 1 for p ∈ {0, 1, 2}. Since the case p = 3 requires additional manipulations, we postpone the analysis of this case to the next Section 5. Note, however, that some of the technical lemmas contained in this present section will also be used in Section 5 to treat the case p = 3. Our main theorem in this section is the following.

Theorem 2. Let g ∈ H 3 (Ω + ), n 2 ∈ C 3,1 (Ω + ) and Γ of class C 4,1 . There are constants ε 0 > 0 and C(p) > 0, independent of ε ∈ (0, ε 0 ], such that u ε -u ε,p H 1 (Ω+) ≤ C(p)ε p+1 , for p ∈ {0, 1, 2} and 0 < ε ≤ ε 0 . (20)
The proof of this theorem requires some preparation. We recall the definition of ũε,3 in (16), the cut off function χ and the coordinates x = (x Γ , ν) in Ω -, both introduced in Section 3, and define ũε,p

χ (x) =      p ℓ=0 ε ℓ u ℓ + (x), x ∈ Ω + , χ(x) p ℓ=0 ε ℓ u ℓ -(x Γ , q/ε), x ∈ Ω -, p ∈ {0, 1, 2, 3}. (21)
Note that ũε,p χ precisely captures the expansions in ( 6) and [START_REF] Duruflé | Higher order generalized impedance boundary conditions in electromagnetic scattering problems[END_REF]. The proof of Theorem 2 is broken into two main steps, since we are not able to show (20) directly, but rather introduce ũε,3 as an intermediate term and prove that

u ε -u ε,p H 1 (Ω+) ≤ u ε -ũε,p H 1 (Ω+) + ũε,p -u ε,p H 1 (Ω+) ≤ Cε p+1 . (22)
To treat the two differences which we need to estimate in the last equation, let us first introduce a sequel of technical lemmas.

Lemma 3. There is a constant C > 0 such that

u 2 L 2 (Γ) ≤ C ∇u L 2 (Ω±) u L 2 (Ω±) + u 2 L 2 (Ω±)
for all u ∈ H 1 (Ω ± ).

Proof. Estimates of this kind are well known from earlier works on generalized impedance boundary conditions, however, since Ω ± is unbounded we briefly sketch the proof for the domain

Ω + . For u ∈ C ∞ (Ω + ) ∩ H 1 (Ω + ), the fundamental theorem of calculus implies |u(x, h)| 2 -|u(x, f (x))| 2 = 2 h f (x)
u(x, x m ) ∂u(x, s) ∂x m ds.

INRIA

If we define the segment S x := {y = (x, s), f (x) < s < h}, then

|u(x, f (x))| 2 ≤ 2 u(x, •) L 2 (Sx) ∂u(x, •) ∂x m L 2 (Sx) + |u(x, h)| 2 .
Integration with respect to h yields

|h -f (x)| |u(x, f (x))| 2 ≤ 2|h -f (x)| u(x, •) L 2 (Sx) ∂u(x, •) ∂x m L 2 (Sx) + u(x, h) 2 L 2 (Sx) .
We integrate with respect to x,

R m-1 |u(x, f (x))| 2 dx ≤ 2 u L 2 (Ω+) ∂u ∂x m L 2 (Ω+) + sup x∈R m-1 |h -f (x)| -1 u 2 L 2 (Ω+) ,
which gives the claim of the lemma, since

u 2 L 2 (Γ) ≤ C R m-1 |u(x, f (x))| 2 dx and since C ∞ (Ω + )∩ H 1 (Ω + ) is dense in H 1 (Ω + ). Lemma 4. Assume that v ε ∈ H 1 (Ω) satisfies ∆v ε + k 2 n 2 v ε = 0 in Ω, ∂v ε ∂ν = T + n 2 + (v ε ) on Γ a for ε ∈ (0, ε 0 ]
and the a-priori estimate

Ω |∇v ε | 2 -k 2 n 2 |v ε | 2 dx - Γa v ε T + n 2 + (v ε ) ds ≤ C ε s+1/2 v ε L 2 (Γ) + ε s v ε L 2 (Ω-) (23)
for C, s > 0 independent of ε and v ε . Then

v ε H 1 (Ω) ≤ Cε s+1 , v ε L 2 (Ω-) ≤ Cε s+2 , and v ε L 2 (Γ) ≤ Cε s+3/2 for ε ∈ (0, ε 0 ]. Proof. Suppose that λ ε = v ε L 2 (Ω) /ε s+1 is unbounded as ε → 0. We set w ε = v ε / v ε L 2 (Ω)
and note that (23) implies

Ω |∇w ε | 2 -k 2 n 2 |w ε | 2 dx - Γa w ε T + n 2 + (w ε ) ds ≤ C λ ε ε -1/2 w ε L 2 (Γ) + ε -1 w ε L 2 (Ω-) . (24)
As C/λ ε is bounded as ε → 0, the last inequality implies an estimate for the absolute value of the imaginary part of the left hand side. Moreover, as -k 2 Im(n 2 ) ≤ 0 and -

Im Γa w ε T + n 2 + (w ε ) ds ≤ 0 as well, we obtain k 2 Ω- |w ε | 2 dx ≤ C λ ε ε 3/2 w ε L 2 (Γ) + ε w ε L 2 (Ω-) .
By Lemma 3,

w ε 3/2 L 2 (Ω-) ≤ Cε 3/2 ∇w ε 1/2 L 2 (Ω-) + w ε 1/2 L 2 (Ω-) w ε 1/2 L 2 (Ω-) + Cε w ε L 2 (Ω-) , RR n°6841
and

w ε L 2 (Ω-) ≤ Cε 3/2 ∇w ε 1/2 L 2 (Ω-) + Cε w ε 1/2 L 2 (Ω-) ≤ Cε 3/2 ∇w ε 1/2 L 2 (Ω-) + Cε 3/2 + 1 2 w ε 3/2 L 2 (Ω-) ,
where we used Young's inequality in the last step. We conclude that

w ε 3/2 L 2 (Ω-) ≤ Cε 3/2 1 + ∇w ε 1/2 L 2 (Ω-) . (25)
Since -Re Γa w ε T + n 2 + (w ε ) ds > 0, the corresponding estimate for the real part of (24) yields

Ω |∇w ε | 2 dx ≤ k 2 Ω Re(n 2 )|w ε | 2 dx + C λ ε ε -1/2 w ε L 2 (Γ) + ε -1 w ε L 2 (Ω-) ≤ C + C ε -1/2 w ε L 2 (Γ) + ε -1 w ε L 2 (Ω-) . ( 26
)
The term ε -1/2 w ε L 2 (Γ) can be estimated by Lemma 3,

ε -1/2 w ε L 2 (Γ) ≤ 2ε -1/2 w ε L 2 (Ω-) + 2ε -1/2 w ε 1/2 L 2 (Ω-) ∇w ε 1/2 L 2 (Ω-) ,
which we plug into (26) to obtain

Ω |∇w ε | 2 dx ≤ C + Cε -1 w ε L 2 (Ω-) 1 + ∇u ε 1/2 L 2 (Ω-) , for ε ∈ (ε, ε 0 ]. ( 27
)
By means of (25),

ε -1 w ε L 2 (Ω-) ≤ C 1 + ∇u ε 1/3 L 2 (Ω-) . (28)
Hence, using the latter bound in (27) we obtain

∇w ε 2 L 2 (Ω) ≤ C(1+ ∇w ε 2/3 L 2 (Ω-) ), which shows that ∇w ε 2 L 2 (Ω) is uniformly bounded in ε ∈ (0, ε 0 ]. From (25) we conclude that ε -1 w ε L 2 (Ω-)
is bounded as well. Lemma 3 implies that ε -1 w ε 2 L 2 (Γ) is bounded for ε ∈ (0, ε 0 ]. Then, however, we conclude from Lemma 12 that ε -1 w ε 2 L 2 (Ω+) is bounded for ε ∈ (0, ε 0 ]. Using (28) we see that ε -1 w ε 2 L 2 (Ω-) is bounded as well. However, the conclusion that ε -1 w ε 2 L 2 (Ω) is bounded for ε ∈ (0, ε 0 ] is a contradiction since by construction w ε 2 L 2 (Ω) = 1. Hence, the bound v ε L 2 (Ω) ≤ Cε s+1 holds for ε ∈ (0, ε 0 ]. The three estimates stated in the lemma now follow from the above inequalities precisely as in the proof of [START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF]Lemma 5.3]. Now we proceed with the estimate of the first difference u εũε,p H 1 (Ω+) appearing in (22).

Theorem 5. Assume that g ∈ H 3 (Ω + ) and n 2 ∈ C 4,1 (Ω + ). Then there is a constant C(p) independent of ε such that

u ε -ũε,p H 1 (Ω+) ≤ C(p)ε p+1
for ε ∈ (0, ε 0 ] and p = {0, 1, 2, 3}.

Proof. We will first prove an estimate for the difference u εũε,p χ H 1 (Ω) , where ũε,p χ has been defined in (21). Since ũε,p χ = ũε,p in Ω + , the estimate stated in the theorem will follow. In the proof of Lemma 5.1 in [START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF] the authors show that the error e ε p := u εũε,p χ satisfies the following transmission problem,

∆e ε p + k 2 n 2 e ε p = 0 in Ω + , ∂e ε p ∂ν = T + n 2 + (e ε p ) on Γ a , [u] Γ = 0, ∂u ∂ν Γ = ε p ∂u p + ∂ν Γ , ∆e ε p + k 2 n 2 e ε p = r p,ε in Ω -, (29)
INRIA where r ε,p (x) = -ε p-1 χ 8 ℓ=1 ℓ-1 j=0 ε p A ℓ-j-1 u k+j+1-ℓ - (x Γ , q/ε) + 2∇χ • p ℓ=0 ε ℓ ∇u p -+ ∆χ p ℓ=0 ε ℓ ∇u p -, for x = x Γ + qν ∈ Ω -.
The computations leading to the form of r ε,p are literally the same as in [START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF]. Also, the proof that g ε,p L 2 (Ω-) ≤ Cε p-1/2 can be achieved as in that reference. The variational formulation for e ε p in (29), obtained with the help of Green's first identity, is

Ω ∇e ε p • ∇v -k 2 n 2 e ε p v dx - Γa v T + n 2 + (e ε p ) ds = ε p Γ ∂u p + ∂ν v ds - Ω- g ε,p v ds.
for all v ∈ H 1 (Ω). Plugging in v = e ε p we find that

Ω ∇e ε p 2 -k 2 n 2 |e ε p | 2 dx - Γa v T + n 2 + (e ε p ) ds ≤ ε p ∂u p + ∂ν L 2 (Γ) e ε p L 2 (Γ) + g ε,p L 2 (Ω-) e ε p L 2 (Ω+) ≤ Cε p e ε p L 2 (Γ) + Cε p-1/2 e ε p L 2 (Ω+)
due to Lemma 1. From this estimate, we conclude by Lemma 4 that e ε p H 1 (Ω) ≤ Cε p+1 , e ε p L 2 (Ω-) ≤ Cε p+2 and e ε p L 2 (Γ) ≤ Cε p+3/2 , for p ∈ {0, 1, 2, 3, 4}. Therefore we can finally estimate

u ε -ũε,p H 1 (Ω+) = u ε -ũε,p+1 + ε p+1 u p+1 + H 1 (Ω+) = u ε -ũε,p+1 χ + ε p+1 u p+1 + H 1 (Ω+) ≤ u ε -ũε,p+1 χ H 1 (Ω+) + ε p+1 u p+1 + H 1 (Ω+) ≤ Cε p+3/2 + Cε p , p = {0, 1, 2, 3}.
For the latter bound we exploited the bounds given in Lemma 1 for u ℓ + H 1 (Ω+) , ℓ ≤ 4, valid due to our assumption that g ∈ H 3 (Ω + ) and n 2 ∈ C 4,1 (Ω + ).

The next theorem proves existence and uniqueness as well as stability of rough surface scattering problems with (ordinary) impedance boundary conditions on Γ. With this result, we will easily be able to prove the error estimate given in Theorem 2. Proposition 6. Let h ∈ H 1/2 (Γ) and µ ∈ C 1 (Γ) such that -Im(µ) ≥ c 0 |µ| 2 > 0 and Re(µ) ≥ c 0 > 0. The impedance problem

∆u + k 2 n 2 u = g in Ω + , ∂u ∂ν = T + n 2 + (u) on Γ a , u + µ ∂u ∂ν = h on Γ, has a unique variational solution u ∈ H 1 (Ω + ), that is, u solves Ω+ ∇u • ∇v -k 2 n 2 uv dx - Γa v T + n 2 + (u) ds + Γ 1 µ uv ds = - Ω+ gv dx + Γ 1 µ hv ds (30)
for all v ∈ H 1 (Ω + ). Additionally, the bound

u H 1 (Ω+) ≤ C 1 + k 2 n 2 L ∞ (Ω+) |µ| 2 / Im(µ) L ∞ (Γ) 1/µ L ∞ (Γ) h H 1/2 (Γ) + µ C 1 (Γ) g L 2 (Ω+)
holds with C independent of µ.

Actually, existence of a solution to the above problem can also be shown under the less restrictive assumption that h ∈ L 2 (Γ) and µ ∈ L ∞ (Γ) with Im(µ) > c 0 > 0; then, however, we did not arrive at a constant in the stability estimate which is homogeneous in µ.

Proof. Our proof is based on the boundedness of the L 2 solution operator shown in Lemma 12. In the first part of the proof, we assume that g = 0. An application of Green's identity (or, plugging in the impedance boundary condition into the variational formulation (30) with v = u)

Ω+ |∇u| 2 -k 2 n 2 |u| 2 dx - Γa u T + n 2 + (u) ds + Γ µ ∂u ∂ν 2 u ds = Γ h ∂u ∂ν ds. ( 31
)
Take the imaginary part of this equation and note that Im Γa u T + n 2 + (u) ds ≥ 0 as well as Im(µ) ≤ 0 to obtain the following estimate,

Γ Im(µ) ∂u ∂ν 2 u ds ≤ h L 2 (Γ) ∂u ∂ν L 2 (Γ) , hence, if Im(µ) ≥ c 0 > 0, ∂u/∂ν L 2 (Γ) ≤ 1/ Im(µ) L ∞ (Γ) h L 2 (Γ) .
Further, taking the imaginary part of (30),

Im Γa u T + n 2 + (u) ds -Im Γ 1 µ uu ds = Im Γ 1 µ hu ds. Since Im Γa u T + n 2 + (u) ds ≥ 0, - Γ Im(µ -1 )|u| 2 ds = Γ Im(µ) |µ| 2 |u| 2 ds ≤ 1/µ L ∞ (Γ) h L 2 (Γ) u L 2 (Γ) and u 2 L 2 (Γ) ≤ |µ| 2 / Im(µ) L ∞ (Γ) 1/µ L ∞ (Γ) h L 2 (Γ) u L 2 (Γ) . (32)
Lemma 12 implies

u 2 L 2 (Ω+) ≤ C |µ| 2 / Im(µ) L ∞ (Γ) 1/µ L ∞ (Γ) h L 2 (Γ) u L 2 (Γ) .
Taking the real part of (31) and noting that -Re Γa u T +

n 2 + (u) ds ≥ 0 yields an H 1 (Ω + ) a-priori estimate for u, u 2 H 1 (Ω+) + Γ Re(µ) ∂u ∂ν 2 u ds ≤ 1 + k 2 n 2 L ∞ (Ω+) u 2 L 2 (Ω+) + h H 1/2 (Γ) ∂u ∂ν 2 H -1/2 (Γ) ≤ C 1 + k 2 n 2 L ∞ (Ω+) |µ| 2 / Im(µ) L ∞ (Γ) 1/µ L ∞ (Γ) C(k,n 2 ,µ) h L 2 (Γ) u H 1 (Ω+) + C h H 1/2 (Γ) u 2 H 1 (Ω+)
and implies existence and uniqueness of solution due to [9, Theorem 2.15], as well as the bound on u given in the theorem. Of course, this is up to now only valid for g = 0.

For g ∈ L 2 (Ω + ) arbitrary, let w ∈ H 1 (Ω + ) be the variational solution of ∆w + k 2 n 2 w = g in Ω + subject to a Dirichlet boundary condition w = 0 on Γ and ∂w/∂ν = T + n 2 + (w) on Γ a . From Theorem 11 we note that w H 1 (Ω+) + ∂u/∂ν L 2 (Γ) ≤ C g L 2 (Ω+) . Due to regularity result INRIA contained in the same lemma, we even have w

H 2 (Ω+) ≤ C g L 2 (Ω+) . The difference v = u -w solves ∆v + k 2 n 2 v = 0 in Ω + , v + µ(∂v/∂ν) = h -µ(∂w)/(∂ν) and ∂v/∂ν = T + n 2 + (w) on Γ a , again in a variational sense. Therefore v H 1 (Ω+) ≤ C 1 + C(k, n 2 , µ) h -µ(∂w/∂ν) H 1/2 (Γ) ≤ C 1 + C(k, n 2 , µ) h H 1/2 (Γ) + µ H 1 (Γ) u H 2 (Ω+) ≤ C 1 + C(k, n 2 , µ) h H 1/2 (Γ) + µ C 1 (Γ) g L 2 (Ω+) .
An application of the triangle inequality to v = u -w shows a H 1 (Ω + ) a-priori bound for u and finishes the proof, again by [START_REF] Ihlenburg | Finite element analysis of acoustic scattering[END_REF]Theorem 2.15].

Corollary 7. (a)

The solution u ε,0 of (19) exits in H 1 (Ω + ), it is independent of ε > 0 and satisfies the bound

u ε,0 H 1 (Ω+) ≤ C g L 2 (Ω) . (b) The solution u ε,2 of (19) exists in H 1 (Ω + ) for all ε > 0 and satisfies u ε,1 H 1 (Ω+) ≤ C g L 2 (Ω) with C independent of ε ∈ (0, ε 0 ] for arbitrary ε 0 .
(c) There is ε 0 > 0 such that the solution u ε,2 of (19) exists in H 1 (Ω + ) for all ε ∈ (0, ε 0 ] and satisfies the bound

u ε,1 H 1 (Ω+) ≤ C g L 2 (Ω) with C independent of ε ∈ (0, ε 0 ].
Proof. (a) This follows from Theorem 11.

(b) The impedance boundary condition for p = 1 is u ε,1 + (ε/α)(∂u ε,1 /∂ν) = 0 on Γ. Therefore, Theorem 6 holds with |q| 2 / Im(q) = √ 2ε and

|µ| 2 / Im(µ) L ∞ (Γ) 1/µ L ∞ (Γ) = √ 2. (c) For p = 2, the impedance boundary condition in (19) is u ε,2 + (ε/α -iHε 2 )(∂u ε,2 /∂ν) = 0 on Γ. Choose ε 0 such that |Hε| < 1/(2 √ 2) for ε ∈ (0, ε 0 ]. Then, Theorem 6 holds with |q| 2 / Im(q) ≤ √ 2ε, |µ| 2 / Im(µ) L ∞ (Γ) 1/µ L ∞ (Γ) ≤ 2 √ 2 and µ C 1 (Γ) ≤ 2ε + ε 2 ∇H L ∞ (Γ)
. Note here that we assumed Γ to be of class C 2 , such that ∇H L ∞ (Γ) is well defined. Now we finish this section with the proof of the main theorem, Theorem 2, showing the convergence u ε,p → u ε in H 1 (Ω + ) of optimal order p + 1.

Proof of Theorem 2. In view of (22) we only need to show that ũε,p -u ε,p H 1 (Ω+) ≤ Cε p+1 . By (17), the difference e ε,p = ũε,p -u ε,p solves ∆e ε,p + k 2 n 2 e ε,p = 0 in Ω + , ∂e ε,p ∂ν = T + n 2 + (e ε,p ) on Γ a , and e ε,p + D ε,p ∂e ε,p ∂ν = ε p+1 r ε,p on Γ.

where r ε,p is given by (18). By Lemma 1 we observe that r ε,p H 1/2 (Ω) ≤ C(p) for p ∈ {0, 1, 2, 3}. Now, Theorem 7 yields the claim of Theorem 2.

GIBC of Order Three

Let us now start to investigate the Neumann-to-Dirichlet impedance boundary condition of order three, with the aim to prove an analogous convergence result as in Theorem 2 for p = 3. As we saw in Section 3, formal expansion yields the following candidate for an Neumann-to-Dirichlet impedance condition of order three,

u ε,3 + ε √ 2 2 (1 + i) - √ 2iεH -(1 -i) ε 2 2 (3H 2 -G + k 2 + ∆ Γ ) ∂u ε,3 ∂ν = 0 on Γ. ( 33
)
The difficulty with this boundary condition is that the operator applied to the normal derivative is a tangential differential operator of order 2. In contrast, the exact Neumann-to-Dirichlet operator is a pseudo differential operator of order -1, mapping H -1/2 (Γ) continuously into H 1/2 (Γ). Also, we remark that for the boundary condition in (33), the term Im Γ u ε,3 (∂u ε,3 /∂ν) ds is indefinite, a fact which would cause trouble in our later analysis. For these two reasons we modify the above condition in the following way. First considering all real terms of the boundary operator in (33), we formally compute that

ε √ 2 2 1 - ε 2 2 (3H 2 -G + k 2 + ∆ Γ ) = ε √ 2 2 1 - ε 2 2 (3H 2 -G + k 2 ) - ε 2 2 (1 -ε 2 ∆ Γ ) -1 ∆ Γ ) + O(ε 4 ), since ∆ Γ -(I -ε 2 ∆ Γ ) -1∆ Γ = -ε 2 ∆ 2 Γ (I -ε 2 ∆ Γ ) -1 .
We note that 1 -ε 2 ∆ Γ is bounded and coercive on H 1 (Γ), hence invertible due to Lax-Milgram's theorem, and that the term on the right constitutes a tangential differential operator on Γ of order zero. Next, we use a Padé approximation to change the sign of the surface Laplacian appearing in the complex terms of (33). Using that (1 -

ε 2 ∆ Γ )(1 + ε 2 ∆ Γ ) = 1 -ε 4 ∆ 2 Γ , we find ε √ 2 2 1 -ε √ 2H + ε 2 2 (3H 2 -G + k 2 + ∆ Γ ) = ε √ 2 2 1 2 -ε √ 2H + ε 2 2 (3H 2 -G + k 2 ) + 1 2 + ε 2 2 ∆ Γ = ε √ 2 2 1 2 -ε √ 2H + ε 2 2 (3H 2 -G + k 2 ) + 1 2 1 -ε 2 ∆ Γ -1 + O(ε 4 ).
Therefore we define a modified boundary operator Dε,3 for a generalized impedance boundary condition of order three as follows,

Dε,3 φ = ε √ 2 2 1 - ε 2 2 (3H 2 -G + k 2 ) - ε 2 2 (1 -ε 2 ∆ Γ ) -1 ∆ Γ ) φ + iε √ 2 2 1 2 -ε √ 2H + ε 2 2 (3H 2 -G + k 2 ) + 1 2 I -ε 2 ∆ Γ -1 φ.
Using this operator, we define u ε,3 via the following boundary value problem

∆u ε,3 + k 2 n 2 u ε,3 = g in Ω + , ∂u ε,3 ∂ν = T + n 2 + (u ε,3 ) on Γ a , u ε,3 + Dε,3 (∂u ε,3 /∂ν) = 0 on Γ. ( 34
)
Three important properties of Dε,3 are collected in the next lemma, namely boundedness and coercivity on L 2 (Γ), which results of course in invertibility of Dε,3 on L 2 (Γ).

Lemma 8. There are constants ε 0 > 0 and C 1,2,3,4 independent of ε > 0 such that

Dε,3 φ L 2 (Γ) ≤ C 1 ε φ L 2 (Γ) , Re Dε,3 φ, φ ≥ C 2 ε φ 2 L 2 (Γ) , and
Im Dε,3 φ, φ ≥ C 3 ε φ 2 L 2 (Γ) for all φ ∈ L 2 (Γ) and ε ∈ (0, ε 0 ]. ( 35
)

INRIA

Proof. We start with some preparations. First, the identity (I -

ε 2 ∆ Γ )(1 -ε 2 ∆ Γ ) -1 φ = φ for φ ∈ L 2 (Γ) implies by one partial integration that Γ ε 2 ∇ Γ ((1 -ε 2 ∆ Γ ) -1 φ) • ∇ Γ ψ + ((1 -ε 2 ∆ Γ ) -1 φ)ψ ds = Γ (1 -ε 2 ∆ Γ ) -1 φψ ds (36)
for all ψ ∈ H 1 (Γ). Setting ψ = (1 -ε 2 ∆ Γ ) -1 φ and taking the complex conjugate of (36) we obtain

(I -ε 2 ∆ Γ ) -1 φ, φ = (I -ε 2 ∆ Γ ) -1 φ 2 L 2 (Γ) + ε 2 ∇ Γ (I -ε 2 ∆ Γ ) -1 φ 2 L 2 (Γ) (37)
as well as ψ = ∆ Γ (1 -ε 2 ∆ Γ ) -1 results by further integrations by parts in

-∆ Γ (I -ε 2 ∆ Γ ) -1 φ, φ = ∇ Γ (I -ε 2 ∆ Γ ) -1 φ 2 L 2 (Γ) + ε 2 ∆ Γ (I -ε 2 ∆ Γ ) -1 φ 2 L 2 (Γ) . (38)
The Cauchy-Schwarz inequality applied to (36) shows that (I -

ε 2 ∆ Γ ) -1 φ 2 L 2 (Γ) ≤ φ 2 L 2 (Γ) . Therefore (37) implies ε 2 ∇ Γ (I -ε 2 ∆ Γ ) -1 φ 2 L 2 (Γ) ≤ φ 2 L 2 (Γ) . (39)
The identity

-ε 2 ∆ Γ (I -ε 2 ∆ Γ ) -1 = I -(I -ε 2 ∆ Γ ) -1 and (37) finally shows that -ε 2 ∆ Γ (I - ε 2 ∆ Γ ) -1 φ, φ ≤ φ 2 L 2 (Γ)
as well, and hence (38

) yields ε 4 ∆ Γ (I -ε 2 ∆ Γ ) -1 φ 2 L 2 (Γ) ≤ φ 2 L 2 (Γ) .
A brief look at the definition of Dε,3 shows that we can write

Dε,3 = ε a ε - √ 2ε 2 /4 (1 -ε 2 ∆ Γ ) -1 ∆ Γ ) + iε b ε + 1 2 I -ε 2 ∆ Γ -1
and that we observe that there is ε 0 > 0 such that the coefficients a ε and b ε on Γ are bounded between two constants 0

< c 0 < C uniformly in ε ∈ (0, ε 0 ]. Consequently, Dε,3 φ, ψ = ε a ε φ, ψ - √ 2ε 3 4 (1 -ε 2 ∆ Γ ) -1 ∆ Γ φ, ψ + iε b ε φ, ψ - iε 2 I + ε 2 ∆ Γ -1 φ, ψ
for all ψ ∈ L 2 (Γ) and for ψ = φ we can estimate real and imaginary part of this expression from below due to our above computations, c 0 φ 2 L 2 (Γ) ≤ Re Dε,3 φ, φ and c 0 φ 2 L 2 (Γ) ≤ Im Dε,3 φ, φ . For an upper bound on L 2 (Γ), we set ψ = Dε,3 φ, which yields by the Cauchy-Schwarz inequality

Dε,3 φ L 2 (Γ) ≤ Cε φ L 2 (Γ) + Cε 3 (I -ε 2 ∆ Γ ) -1 ∆ Γ φ L 2 (Γ) + Cε (I -ε 2 ∆ Γ ) -1 φ L 2 (Γ) ≤ Cε φ L 2 (Γ) .
In the next step, we show that the approximated scattering problem involving the Neumannto-Dirichlet impedance operator Dε,3 is uniquely solvable. Recall that this scattering problem is to find

u ε,3 ∈ H 1 (Ω + ) such that ∆u + k 2 u = g in Ω + , ∂u ∂ν = T + n 2 + (u) on Γ a , u + Dε,3 ∂u ∂ν = 0 on Γ. ( 40
)
Actually, we are going to tackle this problem in more generality, which will be helpful in the later error analysis, and solve for an inhomogeneous boundary condition u + Dε,3 (∂u/∂ν) = h.

Before stating our result on existence and uniqueness of solution, we note that the variational form of the problem then reads

Ω+ ∇u • ∇v -k 2 uv ds - Γa v T + n 2 + (u) ds + Γ v D-1 ε,3 (u) ds = - Ω+ gv dx + Γ v D-1 ε,3 h ds (41)
for all v ∈ H 1 (Ω + ). Note that D-1 ε,3 is indeed boundedly invertible on L 2 (Γ) for ε in some interval (0, ε 0 ] due to its coercivity. Theorem 9. For g ∈ L 2 (Ω + ) and h ∈ H 1/2 (Γ) there is a unique solution of the variational problem (41) which satisfies

u ε,3 H 1 (Ω+) + C 1 √ ε ∂u ε,3 ∂ν L 2 (Γ) ≤ C 2 g L 2 (Ω) + h H 1/2 (Γ)
for C 1,2 independent of ε ∈ (0, ε 0 ].

Proof. As in the proof of Theorem 6 we consider first the case g = 0; the general case will again follow from this special situation. By an integration by parts in Ω + , we find

Ω+ ∇u ε,3 2 -k 2 n 2 |u ε,3 | 2 dx - Γa u ε,3 T + n 2 + (u ε,3 ) ds - Γ u ε,3 ∂u ε,3 ∂ν ds = 0. ( 42
)
Taking the imaginary part of this equation, we obtain

Im Γa u ε,3 T + n 2 + (u ε,3 ) ds + Im Γ u ε,3 ∂u ε,3 ∂ν ds = 0.
Moreover, exploiting the coercivity of Im Dε,3 proven in Lemma 8 Therefore, by Lemma 12 and the fact that Im Γa u ε,3 T +

n 2 + (u ε,3 ) ds ≥ 0, ε ∂u ε,3 ∂ν 2 L 2 (Γ) ≤ h H 1/2 (Γ) ∂u ε,3 ∂ν H -1/2 (Γ) . ( 44
) Due to Lemma 8, D ε,3 is bounded on L 2 (Γ) with bound C 1 √ ε for ε ∈ (ε, ε 0 ].
Exploiting the impedance boundary condition in (41), we obtain

u ε,3 2 L 2 (Γ) ≤ 2 D ε,3 (∂u ε,3 /∂ν) 2 L 2 (Γ) + 2 h 2 L 2 (Γ) ≤ Cε 2 ∂u ε,3 /∂ν 2 L 2 (Γ) + 2 h 2 L 2 (Γ) ≤ Cε h H 1/2 (Γ) ∂u ε,3 ∂ν H -1/2 (Γ) + 2 h 2 L 2 (Γ) ≤ C h H 1/2 (Γ) u ε,3 H 1 (Ω+) + 2 h 2 L 2 (Γ) ,
INRIA for a different constant C. We take the real part of (42) and exploit that -Re Γa u ε,3

≥ cε ∂u ε,3 ∂ν 2 L 2 (Γ) -Re Γ h ∂u ε,3 ∂ν ds.
With the help of Lemma 12, (42) therefore implies

u ε,3 2 H 1 (Ω+) + cε ∂u ε,3 ∂ν 2 L 2 (Γ) ≤ h H 1/2 (Γ) ∂u ε,3 /∂ν H -1/2 (Γ) + (1 + k 2 n 2 ∞ ) u ε,3 2 L 2 (Ω+) ≤ C 1 h H 1/2 (Γ) u ε,3 H 1 (Ω+) + C 2 u ε,3 2 L 2 (Γ) ≤ C1 h H 1/2 (Γ) u ε,3 H 1 (Ω+) + 2C 2 h 2 L 2 (Γ) ,
where C 1 , C 2 and C1 are constants independent of ε. Hence, we conclude that

u ε,3 H 1 (Ω+) + c √ ε ∂u ε,3 ∂ν L 2 (Γ) ≤ C h H 1/2 (Γ) (45)
for some constant C independent of ε.

Consider now the general case 0 = g ∈ L 2 (Ω + ). Again, let w ∈ H 1 (Ω + ) be the variational solution of ∆w + k 2 n 2 w = 0 in Ω + subject to the homogeneous Dirichlet boundary condition u| Γ = 0 and the radiation condition ∂w/∂ν = T + n 2 + (w) on Γ a . Recall that Theorem 11 states

w H 1 (Ω+) + ∂u/∂ν L 2 (Γ) ≤ C g L 2 (Ω+)
and by the regularity result contained in the same theorem we also note the bound w

H 2 (Ω+) ≤ C g L 2 (Ω+) . The difference v = u -w solves ∆v + k 2 n 2 v = 0 in Ω + , v + Dε,3 (∂v/∂ν) = h -Dε,3 (∂w)/(∂ν) and ∂v/∂ν = T + n 2 + (w) on Γ a , again
in a variational sense. The boundedness of Dε,3 on H 1/2 (Γ) shown in Lemma [START_REF] Haddar | Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case[END_REF] implies

v H 1 (Ω+) + c √ ε ∂u ε,3 ∂ν L 2 (Γ) ≤ C(1 + k 2 n 2 ∞ ) h H 1/2 (Γ) + Dε,3 (∂w)/(∂ν) H 1/2 (Γ) ≤ C(1 + k 2 n 2 ∞ ) h H 1/2 (Γ) + C ∂w/∂ν H 1/2 (Γ) ≤ C(1 + k 2 n 2 ∞ ) h H 1/2 (Γ) + C u H 2 (Ω+) We conclude that u H 1 (Ω+) + c √ ε ∂u ε,3 ∂ν L 2 (Γ) ≤ v H 1 (Ω+) + w H 1 (Ω+) ≤ C(1 + k 2 n 2 ∞ ) h H 1/2 (Γ) + C g L 2 (Ω+) ,
which is H 1 (Ω + ) a-priori bound for u and finishes the proof, since this a-priori bound establishes existence and uniqueness of solution, again by [START_REF] Ihlenburg | Finite element analysis of acoustic scattering[END_REF]Chapter 2].

The latter existence and stability results provides the mean to prove an optimal error estimate for u ε -u ε,3 H 1 (Ω+) . Of course, this estimate makes use of several of lemmas contained in the previous Section 4. The error estimate given here requires smoothness assumptions on the data g, n 2 and Γ of our scattering problem. We give those assumptions explicitly, noting that we could of course also simply require that all data be smooth enough.

Theorem 10. Assume that g ∈ H 6 (Ω + ), n 2 ∈ C 6,1 (Ω) and that Γ is of class C 7,1 . Then there are constants ε 0 > 0 and C > 0, independent of ε ∈ (0, ε 0 ], such that

u ε -u ε,3 H 1 (Ω+) ≤ Cε 4 . (46)
Proof. As in (22) we split the quantity we need to estimate in two parts,

u ε -u ε,3 H 1 (Ω+) ≤ u ε -ũε,3 H 1 (Ω+) + ũε,3 -u ε,3 H 1 (Ω+) .
The first part has been already treated in Lemma 5. For the second part, we note that e ε,3 := ũε,3 -u ε,3 solves ∆e

ε,3 + k 2 n 2 e ε,3 = 0 in Ω + , ∂e ε,3 ∂ν = T + n 2 + (e ε,3
) on Γ a , and e ε,3 + Dε,3 ∂e ε,k ∂ν = ε 4 rε,3 on Γ.

The right hand side rε,3 is different from r ε,3 given in (18), since rε,3 relies on the modified boundary operator Dε,3 instead of D ε,3 . More precisely,

ε 4 rε,3 = ε 4 r ε,3 + Dε,3 -D ε,3 ∂ ũε,3 ∂ν .
From Lemma 1 it follows that r ε,3 H 1/2 (Γ) is bounded by some constant independent of ε. However, the difference Dε,3 -D ε,3 is, modulo constants independent of ε, sum of the two terms

ε 2 ∆ Γ -ε 2 (I -ε 2 ∆ Γ ) -1 ∆ Γ = -ε 4 ∆ 2 Γ (I -ε 2 ∆ Γ ) -1 and 1 + ε 2 ∆ Γ -(1 -ε 2 ∆ Γ ) -1 = ε 4 ∆ 2 Γ (1 -ε 2 ∆ Γ ) -1 . Thus, to show that ( Dε,3 -D ε,3 )(∂ ũε,3 /∂ν) H 1/2 (Γ) ≤ C for C independent of ε, we need to prove that ∆ 2 Γ (1 -ε 2 ∆ Γ ) -1 (∂ ũε,3 /∂ν) H 1/2 (Γ)
is bounded independently of ε. Therefore it is sufficient to show that this property holds for ∆ 2 Γ (∂ ũε,3 /∂ν) H 1/2 (Γ) , that is, we need that ∂ ũε,3 /∂ν H 9/2 (Γ) is bounded. However, Lemma 1 shows that such a bound is guaranteed if only g, n 2 and Γ are smooth enough; for our purpose, we need for instance g ∈ H 6 (Ω + ). If all data is smooth enough, an application of Theorem 9 finishes the proof.

Preliminary Numerical Results

This section is dedicated to the presentation of some preliminary numerical results that test the accuracy of derived GIBCs. These tests are restricted to the 2-D case. We opted for the use of finite element methods to compute the numerical solution using FreeFem++ library (http://www.freefem.org/ff++). The major difficulty linked to the use of volumic methods is indeed the truncation of the computational domain. One possible approach is the use of Dirichlet-to-Neumann maps introduced in the second section of this paper. The second possibility, adopted here, is the use of so-called perfectly matched absorbing layers (PML) (see Fig. 2). We refer to [START_REF] Chandler-Wilde | The pml for rough surface scattering[END_REF] for description of this method in the context of rough surface scattering problems.

The numerical results presented in the following intends to only give first numerical validations and also hints on possible numerical difficulties that need to be addressed in future works. The two examples, shown in Fig. 3 and 4 correspond to a computational domain which is 100 λ large and 5 λ thick, where λ = 2π/k (= 1 in our cases). The PML width is λ and we used an (optimized) linear profile for the absorption coefficient going from 1 at the inner boundary of the PML to σ 0 = 5.75 + 3i at the outer boundary (see [START_REF] Chandler-Wilde | The pml for rough surface scattering[END_REF]). The numerical examples are computed using P 2 finite elements. A reference solution is computed using a fine mesh of the two layered medium.

We observe in both experiments how the use of higher order GIBC significantly improves the accuracy of the approximate model. For all GIBCs, the error decreases with respect to ε. However, we failed in obtaining confirmation of convergence rates as predicted by our theory. We believe that this is mainly due to the numerical error induced by lateral bounds of the computational domain (the PML method can only be justified if an infinite horizontal layer is used above the rough surface). We observed that the residual error due to this truncation cannot be made smaller than 1%, which would also explain why the three curves (for GIBCs) meet for small ε. We therefore think that the design of efficient method to bound the computational domain constitutes one important step before going deeper in the numerical validation of these approximate boundary conditions. Addressing this issue is far beyond the scope of the present work.

We end us this discussion by noticing that when comparing the two examples, one observes that better results are obtained for rough surfaces with less sharp variations, which is somehow a predictable behavior since more regularity is needed for higher order GIBCs.

A Auxiliary Results on Rough Surface Scattering

In this section, we provide some results on existence and uniqueness for Dirichlet scattering problems in Ω + . The following Rellich identity, which is valid for any solution v in H

u H 1 (Ω) + ∂u/∂ν L 2 (Γ) ≤ C g L 2 (Ω) , (48)
and elliptic regularity implies u H 2 (Ω) ≤ C g L 2 (Ω) . In the general case,

u H 1 (Ω) ≤ C g L 2 (Ω) + h H 1/2 (Γ) .
Generally

Figure 1 :

 1 Figure 1: The geometry for the rough layer problem. The domain Ω = Ω + ∪ Γ ∪ Ω -lies in between the two planes Γ a = {x m = a} and Γ -a = {x m = -a}. In Ω, the refractive index n 2 varies; while n 2 is Lipschitz continuous in Ω ± , the index may jump across Γ. The unit normal ν points out of Ω and on Γ we choose ν to point downwards. The domain Ω R = {x ∈ Ω, |x| < R} is obtained from Ω by cut off in the lateral variables x

  1 ), allows to explicitly compute a Dirichlet-to-Neumann operator, mapping Dirichlet values φ on Γ a to the Neumann boundary values of the unique radiating solution u in U + a taking Dirichlet trace values φ on Γ a . Construction of this operator relies on the following representation formula for u ε ,

Figure 2 :

 2 Figure 2: The numerical experiments settings

Figure 3 :Figure 4 :

 34 Figure 3: Left: surfcae profile: f (t) = 3.5 + 0.2 cos(t -3) + 0.4 exp(-(t -3) 2 ). Right: L 2 error versus ε (log-log scale) for GIBCs of order 0, 1, 2 and 3.

2 + 2 ∂v ∂x m 2 + k 2 (x m + a)ν m ∂n 2 ∂x m |v 2 | 2 dx + 2h Γa |∇v| 2 -2 ∂v ∂ν 2 -k 2 n 2 |v| 2 Theorem 11 .+ n 2 +

 22222222112 2 (Ω + ) of the Helmholtz equation ∆v + k 2 n 2 v = g, subject to the radiation condition ∂v/∂ν = T + n (v) on Γ a , irrespectively of the actual boundary condition on Γ, has been shown in[START_REF] Lechleiter | Scattering of acoustic waves from rough layers[END_REF],Ω+ + a) ν m |∇v| 2 -2 Re ∂v ∂x m ∂v ∂ν -ν m k 2 |v| 2 ds -Assume that n 2 ∈ C 0,1 (Ω + ), ∂n 2 /∂x m ≤ 0 in Ω + and that Γ is of class C 2 .Then there is unique variational solution of ∆u + k 2 n 2 u = g in Ω + subject to the Dirichlet boundary condition u = h on Γ, that is, u satisfiesΩ+ ∇u • ∇v -k 2 n 2 uv dx -Γa vT (u) ds = -Ω+ gv dx for all v ∈ H 1 0 (Ω + ) := {v ∈ H 1 (Ω + ), v =0 on Γ} and the boundary condition u| Γ = h in the trace sense. If h = 0, the solution satisfies

2 +∂ν 2 L 2

 222 Ω+gw dx ≤ 2k g L 2 (Ω+) w L 2 (Ω+) .The inequality after the line break follows from taking the imaginary part ofΩ+ |∇w| 2 -k 2 n 2 |w| 2 dx -(x m + a)ν m < 0 and -Re Γa vT + n (v) ds ≥ 0 we obtain that ∂w/(Γ) ≤ C g L 2 (Ω+) w H 1 (Ω+) ≤ C g 2 L 2 (Ω+) .(50)Let us denote Ω + R := Ω + ∩ Ω R . Twice applying Green's first identity yields 0 =Ω + R (∆u + k 2 n 2 u)w dx = Ω + R ∇u • ∇w -k 2 n 2 dsu ∈ H 1 (Ω +) and w ∈ H 2 (Ω + ), both volume integrals converge, as R → ∞. Moreover,

  T +

	0 by the representation of T + n 2 +	in (3), and that						n 2 +	(u ε,3 ) ds ≥
	-Re	Γ	u ε,3 ∂u ε,3 ∂ν	ds = Re	Γ	Dε,3	∂u ε,3 ∂ν	∂u ε,3 ∂ν	ds -Re	Γ	h	∂u ε,3 ∂ν	ds

  , elliptic regularity yields thatu H p+2 (Ω) ≤ C g H p (Ω) + h H p+3/2 (Γ) ,under the smoothness assumption n 2 ∈ C p,1 (Ω + ) and Γ of class C p+1,1 . Due to Theorem 11 there is a unique weak solution u ∈ H1 (Ω + ) of the problem ∆u + k 2 n 2 u = 0 in Ω + , (u) on Γ a , u = h on Γ,for all h ∈ H 1/2 (Γ), and u H 1 (Ω+) ≤ C h H 1/2 (Γ) for C independent of u. Moreover, Theorem 11 states that for arbitrary g ∈ L 2 (Ω + ) there is a unique solution of the problem∆w + k 2 n 2 w = g in Ω + , H 2 (Ω+) ≤ C g L 2 (Ω+) ,again for C independent of u. Moreover, by the Rellich identity (47) and the Dirichlet boundary condition w = 0 on Γ which implies that ∂w/∂x m = e m • (ν∇w) = ν m (∂w/∂ν), + k 2 (x m + a)ν m ∂n 2 ∂x m |w 2 | 2 dx -

								∂u ∂ν	= T + n 2 +
								∂w ∂ν	= T + n 2 +	(u) on Γ a ,	w = 0 on Γ,
	with norm bound w Ω+ 2 ∂w ∂x m 2										Γa	wT + n 2 +	(w) ds -	Γ	(x m + a)ν m	∂w ∂ν	2	ds
	+ 2h	Γa	|∇w| 2 -2	∂w ∂ν	2	-k 2 n 2 |w| 2 ds = -2	Γ	(x m + a) Re	∂w ∂x m	g dx -Re	Ω+	gw dx.
	Due to [3, Lemma 2.2] and since n 2	Γa = n 2	it holds
	2h	Γa	-|∇w| 2 + 2	∂w ∂ν

INRIA

Proof of Lemma 12.

2 + k 2 n 2 + |w| 2 ds ≤ 2k Im

RR n°6841

Proof. The existence statement for h = 0 in H 1 0 (Ω + ) is shown as in [START_REF] Chandler-Wilde | Existence, uniqueness, and variational methods for scattering by unbounded rough surfaces[END_REF], where a Rellich identity similar to (47) is used to prove existence of solution in H 1 (Ω + ). For h = 0, the regularity statement concerning u L 2 (Γ) also stems from the identity (47); in fact, this is shown explicitly in the proof of Lemma 12 below. The existence result for arbitrary h is shown as in the standard proof for inhomogeneous elliptic boundary value problems [START_REF] Mclean | Strongly Elliptic Systems and Boundary Integral Operators[END_REF]Theorem 4.10], noting that (48) even implies u H 1 (Ω) ≤ C g (H 1 0 (Ω+)) * , where (H 1 0 (Ω + )) * denotes the space of bounded linear functionals on (H 1 0 (Ω + )) * with obvious norm. The regularity statement in the general case follows from the corresponding regularity result for bounded domains [START_REF] Mclean | Strongly Elliptic Systems and Boundary Integral Operators[END_REF]Theorem 4.19] and a technique already introduced in [START_REF] Lechleiter | Scattering of acoustic waves from rough layers[END_REF]. For shake of completness we shall hereafter sketch some details of the proof. By abuse of notation we do not distinguish between a solution u ∈ H 1 (Ω + ) of the Helmholtz equation satisfying ∂u/∂ν = T + n 2 + (u) on Γ a and its unique radiating extension

By Q 2 j := j + 2Q we denote an even larger cube containing Q j . From the boundary regularity result [START_REF] Mclean | Strongly Elliptic Systems and Boundary Integral Operators[END_REF]Theorem 4.28] we obtain the following estimate in each cube Q j ,

The constants C j depend on p, the local smoothness of the coefficient n 2 and the local regularity of Γ. Hence, the lemma's assumptions that Γ be of class C p+1,1 and n 2 be C 0,1 imply a uniform bound C for the numbers C j . Now we can exploit that the Q j cover Ω,

The regularity result can then easily derived from (48) and (49) by a simple induction.

Several times in this text, for instance in the proof of Theorem 9, we rely on the proof of the following auxiliary result.

Lemma 12 (L 2 solution operator). For a weak solution u ∈ H 1 (Ω + ) of the Helmholtz equation ∆u + k 2 n 2 u = 0 in Ω + which takes boundary values φ ∈ H 1/2 (Γ) and satisfies the radiation condition ∂u/∂ν = T +

For bounded domains, this is a well known result [START_REF] Mclean | Strongly Elliptic Systems and Boundary Integral Operators[END_REF]Theorem 4.25]. One possibility how to prove Lemma 12 is indeed to mimic the proof of the result for bounded domains, exploiting results from variational theory of rough surface scattering. The other option is to use results from integral equation theory for rough surface scattering problems, namely invertibility in L 2 (Γ) of the boundary integral operator arising in a Brackhage-Werner ansatz for the rough surface Dirichlet scatter problem. This approach is of course restricted to a constant index of refraction in Ω + , and therefore we prefer the first option in the sequel.

The two last terms tend to zero as R → ∞ and the limit for the remaining two terms on the right exist. Hence We set g = u, take the real part of the latter equation and obtain with the help of (50)

which is precisely the claim of the lemma.


Stéphanie Roux

Mathilde Courel

Inès Birlouez-Aragon

Francesco Municino

Mario Massa

Roux Stéphanie

Courel Mathilde

Birlouez-Aragon Inès

Comparative thermal impact of two UHT technologies, continuous ohmic heating and direct steam injection, on the nutritional properties of liquid infant formula

Keywords: Maillard reaction, thermal degradation, kinetics, thermal treatment Highlights

à 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.

Comparative thermal impact of two UHT technologies, continuous ohmic heating and direct steam injection, on the nutritional properties of liquid infant formula

Introduction

Milk based infant formulas (IF) which composition is close to mature human milk, are designed to cover the nutritional needs of healthy infants. Product stability and microbiological safety are usually achieved by atomization or heat sterilization. However, severe heat treatments have deleterious effects on the sensory, biophysical and nutritional properties of milk products [START_REF] Finot | The extent of the Maillard reaction during the processing of milk[END_REF]. Major changes are detected in thermally treated milks like protein denaturation and aggregation [START_REF] Oldfield | Kinetics of Denaturation and Aggregation of Whey Proteins in Skim Milk Heated in an Ultra-high Temperature (UHT) Pilot Plant[END_REF][START_REF] Oldfield | Heat-induced interactions of βlactoglobulin and α-lactalbumin with the casein micelle in pH-adjusted skim milk[END_REF], lipid-protein and protein-protein interactions [START_REF] Fenaille | Modifications of milk constituents during processing: a preliminary benchmarking study[END_REF], sugar isomerization [START_REF] Berg | Degradation of lactose during heating of milk. 1. Reaction pathways[END_REF] and a wide range of chemical reactions, the Maillard reaction being the main involved [START_REF] Hiller | Functional properties of milk proteins as affected by Maillard reaction induced oligomerisation[END_REF][START_REF] Pellegrino | Coupling of lactulose and furosine indices for quality evaluation of sterilized milk[END_REF][START_REF] Van Boekel | Effect of heating on Maillard reactions in milk[END_REF]. It is a complex series of reactions in which, in milk, a condensation reaction first occurs between the carbonyl group of lactose and the -amino group of lysine. The new compounds produced by Maillard reaction are likely to affect the quality, nutritional value and safety of the product [START_REF] Birlouez-Aragon | Milk fortification and processing -Nutritional impact[END_REF].

Two classical modes of high temperature short time (HTST) heating -indirect and direct UHT -are commonly used for sterilization of milk and milk products. In the indirect mode, the milk is separated from the heating fluid (steam or warm water under pressure) by a wall while in the second case it is mixed with the warming fluid (steam). Some of the steam is condensed, giving up its latent heat of vaporization to the product and giving a much more rapid rate of heating than is available with any indirect system. After holding at the sterilizing temperature, the product is expanded in a cooling vessel, where the pressure is at a level below atmospheric. It immediately boils and gives up water in the form of water vapor which is removed from the vessel and condensed. Steam injection (SI) has been observed to be the best technology to limit thermal damage of milks [START_REF] Birlouez-Aragon | A rapid fluorimetric method to estimate the heat treatment of liquid milk[END_REF][START_REF] Van Asselt | Extreme high-temperature treatment of milk with respect to plasmin inactivation[END_REF]. Although it is a very efficient heating method, where the product can reach temperatures above 150 °C in less than a second, the main disadvantage of SI relies on increased complexity and costs [START_REF] Burton | Chapter 5. Performance of UHT processing plants, Ultra-high-temperature processing of milk and milk products[END_REF]. Firstly, the water used for vapor production needs to be of "culinary" grade: it fulfils very strict hygiene standards and the use of chemical sterilizing agents that might end in the food product is a critical point. Boilers and steam generation equipment have to be operated in a way to prevent foaming, priming, carryover and excessive entrainment of boiling water into the steam (Canadian Food Inspection Agency, 2014). Secondly, the dilution of the product caused by steam condensation (up to 11 % of the product) needs to be compensated by a vaporization step to recover the exact initial solid content of the milk; the recovered vapor requires expensive devices to be recycled in comparison with an indirect heat exchanger allowing 80 to 90% heat recovery [START_REF] Burton | Chapter 5. Performance of UHT processing plants, Ultra-high-temperature processing of milk and milk products[END_REF]. Thirdly, elevated costs of heat and water are involved in direct UHT: the steam distribution through a network of pipes down to the injector comes along with heat losses and large amounts of water are required for the operating of the condenser. Furthermore, the residence time distribution can be quite dispersed within the steam nozzle causing thermal heterogeneity. Finally, fouling in the downstream area of the SI heater is an important operating problem, threatening the quality of product and affecting plant operation [START_REF] Truong | The Use of a Heat Flux Sensor for In-Line Monitoring of Fouling of Non-Heated Surfaces[END_REF].

Another direct way of warming up a food product is ohmic heating. It is based on a simple mechanismthe Joule effectwhich consists of raising the product temperature by passing an electrical current directly through the food when placed between two electrodes. Ohmic heating (OH) offers the big advantage over conventional indirect heating methods, of being based on purevolume heating: the heat is directly dissipated in the product which behaves as a resistance thanks to the presence of free ions, resulting in fast and uniform heating. Even though it is based on the use of electricity which is less and less appreciated today, the electrical conversion into heat energy is close to 100% [START_REF] Ghnimi | Energy efficiency of an ohmic heating technology by fluid jet in food industrial process[END_REF] and the very low thermal inertia makes possible fast and precise regulation. Heat transfer is a function of the electrical and thermal parameters of the product and equipment which makes the process easy-to-drive [START_REF] Goullieux | Chapter 22 th Ohmic Heating[END_REF][START_REF] Roux | Design of an ohmic heating reactor to study the kinetics of thermal reactions in liquid products[END_REF]. However, fouling on electrode surfaces constitutes one of the main problems of OH [START_REF] Ayadi | Continuous ohmic heating unit under whey protein fouling[END_REF][START_REF] Ayadi | Experimental study of hydrodynamics in a flat ohmic cell -impact on fouling by dairy products[END_REF][START_REF] Fillaudeau | Heat treatment of whole milk by the direct Joule effect -Experimental and numerical approaches to fouling mechanisms[END_REF][START_REF] Stancl | Milk fouling at direct ohmic heating[END_REF]. Besides production shut down for cleaning, fouling induces quality losses of the product and presence of electric arcs. Two big ways have been explored to limit fouling in OH: the one goes through cooling the walls of the ohmic reactor [START_REF] Pain | Continuous Heat Treatment Method and Heating Device for an Electrically Conductive Fluid[END_REF]; and the second is based on the use of a fluid jet subjected to the electric field [START_REF] Ghnimi | Design and performance evaluation of an ohmic heating unit for thermal processing of highly viscous liquids[END_REF]. The first approach was used in this study to limit fouling problems during sterilization of the liquid infant formula.

The question of quantifying the healthiness of a product, which is crucial for infant formulas, and the possibility of optimizing food processing with respect to health aspects were addressed by [START_REF] Van Boekel | Product quality and food processing: how to quantify the healthiness of a product[END_REF] who observed among others, the (anti)mutagenicity of Maillard products in heated milk. Different authors have used various markers of heat damage in milk and infant formulas in order to evaluate the effect of a given technology on the product quality: the combination of furosine and lactulose values was used by Pellegrino et al. to build information about the thermal history of sterilized milk [START_REF] Hiller | Functional properties of milk proteins as affected by Maillard reaction induced oligomerisation[END_REF]. [START_REF] Rufián-Henares | Generation of Furosine and Color in Infant/Enteral Formula-Resembling Systems[END_REF] monitored the extent of Maillard reaction by measuring furosine and color in model infant formula mixed and heated at lab scale. In addition to furosine, hydroxymethylfurfural (HMF) and pyrraline were studied by [START_REF] Contreras-Calderon | Indicators of nonenzymatic browning in the evaluation of heat damage of ingredient proteins used in manufactured infant formulas[END_REF] as indicators of thermal damage together with available lysine as nutritional indicator, to evaluate heat damage to ingredients used in commercial infant formulas. [START_REF] Morales | Ratio of maltose to maltulose and furosine as quality parameters for infant formula[END_REF] proposed to use the ratio of maltose to maltulose and furosine as quality parameters for infant formula. Rapid fluorimetric methods were investigated to assess the quality of heat-treated food products: the FAST index (Fluorescence of Advanced Maillard products and Soluble Tryptophan) was first designed for liquid milk products [START_REF] Birlouez-Aragon | A rapid fluorimetric method to estimate the heat treatment of liquid milk[END_REF] and extended to liquid infant formulas [START_REF] Birlouez-Aragon | Assessment of protein glycation markers in infant formulas[END_REF][START_REF] Damjanovic Desic | The FAST index -A highly sensitive indicator of the heat impact on infant formula model[END_REF]; it is based on the fluorescence of Maillard products observed in the soluble extract of the food and corrected for the protein concentration of the obtained solution. Front face fluorescence finger prints were used as rapid predictor of the nutritional quality of a bigger range of food products, by correlating the fluorescence response to the concentration of specific Maillard products like carboxymethyllysine (CML) in IF [START_REF] Birlouez-Aragon | Fluoresence fingerprints as a rapid detector of the nutritional quality of processed and stored foods[END_REF], HMF in milk [START_REF] Schamberger | Evaluation of front-face fluorescence for assessing thermal processing of milk[END_REF], or the concentration of vitamin C, protein denaturation and accumulation of Maillard products in IF [START_REF] Diez | Potential of front face fluorescence associated to PLS regression to predict nutritional parameters in heat treated infant formula models[END_REF]. Finally, [START_REF] Feinberg | Evaluation of tracers for the authentication of thermal treatments of milks[END_REF] demonstrated that no tracer can universally discriminate pasteurization from high pasteurization, direct UHT, indirect UHT, and sterilization of commercial milks. He recommended a multivariate approach by combining at least five tracers, the most discriminative ones being those which globally measure the structural modifications of the milk protein rather than those which specifically quantify the metabolites of the Maillard reaction.

In this study two types of continuous direct UHT treatments were applied to sterilize a model infant formula. The first objective was to compare the performances of OH to SI which is considered the best technology for UHT treatments of milk and milk products. Both technologies were applied at pilot scale, switching from one to another on the same pilot plant to generate identical time-temperature histories of the sterilized IF. They were compared on the basis of thermal degradation of proteins or vitamin C and production of Maillard products. The second objective was to generate kinetic data under real HTST processing conditions, including heating, holding and cooling phase, by modulating the holding time and temperature.

4/16

Materials and methods

Liquid infant formula

The liquid infant formula was specially designed for these experiments by industrial partners of ICARE project. It was produced at industrial scale using microfiltrated skim milk supplemented with various other ingredients in order to fulfil the nutritional target of the designed product (Table 1). The liquid product was packed into vacuum bags of 1 ton (2 units, one for each trial), transported and stored at 4 °C for 4 days. (g/L) 1026 ± 4 Viscosity (b) (mPa.s) exp(0.633 -0.016×T + 1.78E-5×T 2 ) Specific heat capacity (c) (J/kg/K) 3939 + 0.16×T + 0.0041×T 2 Thermal conductivity (c) (W/m/K) 0.534 + 0.0016×T -6.37E-6×T 2 Electrical conductivity (d) (S/m) 0.1082 + 0.008142×T (a) Measured (n = 192), decrease of about 6% between 20 and 140 °C (second order equation) c (b) According to [START_REF] Mccarthy | Rheology of liquid and semi-solid milk products[END_REF], considered equivalent to milk (0.03-15% fat, 70-135 °C), T in °C (c) According to [START_REF] Singh | Introduction to Food Engineering[END_REF], temperature range of 0-150 °C, T in °C (d) According to [START_REF] Roux | Design of an ohmic heating reactor to study the kinetics of thermal reactions in liquid products[END_REF], temperature range of 20-140 °C, T in °C 2.2. Analytical methods

FAST method

The analytical method used to monitor the Fluorescence of Advanced Maillard products and Soluble Tryptophan (FAST) was the one previously developed and patented [START_REF] Birlouez-Aragon | Method for evaluating the heat treatment to which a proteinic nutrient such as milk is subjected[END_REF][START_REF] Birlouez-Aragon | A new method for discriminating milk heat treatment[END_REF]: first, a soluble extract of the IF was produced by bringing 1 mL sample to pH 4.6 using 50 mL of sodium acetate buffer (0.1 M, pH 4.6). 4 mL of 5/16 supernatant were then filtered through a 0.45 µm nylon filter (VWR, France) and placed in a disposable 4 faces acryl cuve (Sarstedt, France). The fluorescence of two complementary indicators was measured in the same sample, using a spectrofluorimeter (CaryEclipse Varian, France): at excitation/emission wavelengths of 290/340 nm, measurement of tryptophan fluorescence (Ftrp) enabled to calculate the soluble protein content and at 330/420 nm measurement of advanced Maillard products (Famp) were performed in the acid soluble supernatant. The FAST index was then calculated from both fluorescence signals, giving an indication of the nutritional damage undergone during the thermal treatment:

100 index FAST   trp amp F F (1)
Maximum variation coefficients of 3.5% were obtained for Ftrp, Famp and FAST index. Ftrp was calibrated with soluble whey proteins (Prolacta 90 from Lactalis, Laval, France) solutions of varying concentration Csoluble proteins from 0.2 to 20 g/L. The calibration curve given by Eq. ( 2) showed linearity between these two quantification limits:

Ftrp = (223.75 ± 0.89)*Csoluble proteins -(42.983 ± 5.691) (2)
2.2.2. Furosine and carboxymethyllysine analysis Furosine and carboxymethyllysine (CML) were analyzed in triplicate according to the method proposed by [START_REF] Charissou | Evaluation of a gas chromatography/mass spectrometry method for the quantification of carboxymethyllysine in food samples[END_REF]: the samples were defatted with ethanol and dichloromethane before centrifugation. A proteinic ring was formed, transferred to another tube and hydrolyzed with concentrated HCl. After vacuum drying, the mix was further dissolved in water, filtered (nylon, 0.22 μm) and dried again. Then, a derivatization step comprised esterification of the carboxylic functions (methanol under strongly acidic conditions) and acylation of the amine functions (anhydrous trifluoroacetic acid). Finally, CML and furosine were quantified by selected ion monitoring on a gas chromatograph coupled with tandem mass spectrometry (Thermo Electron Corporation, Waltham, USA). In order to relate the results to the protein concentration, protein quantification was obtained using the rapid fluorescamine method [START_REF] Yaylayan | A fluorescamine-based assay for the degree of glycation in bovine serum albumin[END_REF]: 20 μL of filtered hydrolysate was made up to 3 mL by borate buffer (0.2 M; pH 8.5) and mixed with 100 μL of fluorescamine (0.15 g/L of acetone). The fluorescence was measured at excitation/emission wavelengths of 390/475 nm on a spectrofluorimeter (Fluoromax-Spex, Jobin-Yvon, Longjumeau, France).

To save time, the concentration of furosine and CML was also predicted by frontal fluorescence. Frontal fluorescence measurements were made by using a sensor developed for ICARE 1 project. The method used by [START_REF] Rizkallah | Front face fluorescence spectroscopy and multiway analysis for process control and NFC prediction in industrially processed cookies[END_REF] for biscuits was applied to samples of thawed infant formula and the model to predict concentration was constructed by reference to the GC-MS-MS method described above. The values predicted by frontal fluorescence were mixed to those obtained by GC-MS-MS and no difference between the values obtained by either chemistry or spectrometry will be made in the rest of the study. The limits of quantification were estimated around 1 mg/100 g proteins for furosine and 1 µg/mg proteins for CML.

Vitamin C

Vitamin C was analyzed by HPLC with fluorescence detection according to the method developed by [START_REF] Tessier | Validation of a micromethod for determining oxidized and reduced vitamin C in plasma by HPLC-fluorescence[END_REF] and adapted to dairy product by Gliguem and Birlouez-Aragon 6/16 (2005). Relative contents (in percentage) were expressed by dividing the peak area by the peak area of the control sample. The quantification domain is comprised between areas of 3×10 5 and 3.5×10 6 with coefficient of variation of repeatability lower than 5%.

Color

Color measurements were done directly on 5 mL fresh or defrost samples with a colorimeter (Luci 100, Dr Lange, Villeurbanne, France) and expressed in the CIELAB system with L* a* b* coordinates. Samples were illuminated with D65 artificial daylight (10° standard angle). Relative variations of 1.5%, 2.5% and 2.0% respectively, were obtained for L*, a* and b* parameters due to colorimeter uncertainty (n = 10). To better characterize the variations in color, each parameter was expressed in Δ values meaning that the initial value of the parameter was subtracted to the measured ones.

Continuous sterilization equipment

The pilot plant was located in the Stazione Sperimentale per l'Industria delle Conserve Alimentari (SSICA) in Parma, Italy. The complete sterilization device was divided into four sections dedicated to preheating, heating, holding and cooling (Fig. 1). A flexible assembling of the different sections of the pilot plant enabled interchanging the heating section from ohmic heating to steam injection without changing the rest of the device. The infant formula (4 °C) was directly poured out of the vacuum pouches into a feeding tank and sent by a pressure pump to a first plate heat exchanger for preheating up to 50 °C. The product then entered the heating section, either ohmic or steam injection, to be heated up to the target temperature (120-140 °C). The target temperature was maintained in the holding section which was equipped with seven valves enabling sampling the product at different treatment times (Table 2). The thermal treatment finally ended in the cooling section with a second plate heat exchanger fed with cold water (almost 10 °C). The whole system was maintained under pressure (4 bars) by means of a counter pressure valve placed downstream to the process.

The ohmic cell used was patented [START_REF] Pain | Continuous Heat Treatment Method and Heating Device for an Electrically Conductive Fluid[END_REF] and specially designed to limit fouling. It was based on the principle developed by [START_REF] Pain | Continuous heat treatment process for fluid foodstuff passing through duct uses electrodes in duct walls and controlling of surface temperature[END_REF] and improved regarding energy savings. The ohmic cell was constituted by two coaxial tubes. The warmer liquid circulating through the intlet path is cooled by the colder liquid circulating through the outlet path without any need for a secondary cooling fluid (Fig. 2). The product was continuously heated by electrodes made of titanium and placed along the longitudinal flowing axis.

The steam injection equipment consisted in a direct inlet of steam into the continuous flow of product obtained by connecting an additional pipe to the principal one. The steam flow was manually set to ensure that the product reached the target temperature.

No PID was used to stabilize the temperature: for OH, the temperature of the product was adjusted by delivering the adequate power with a manual setting and for SI by manually increasing or decreasing the injected steam flow rate. Once the appropriate settings were found, a few minutes were spent to confirm the temperature stability before samples could be taken. The temperature was monitored with sensors of maximum +/-1°C accuracy and remained stable all along the experiments at each targeted temperature.

8/16

Sampling was made by simply opening the valve without any specific pressure control, inducing partial vaporization of the product because of pressure drop. Knowing the flow of liquid product -150 L/hand the length and diameter of the tubes, seven treatment times could be calculated for the ohmic heating device from 1.6 s to 154.3 s and for the steam injection one from 0.4 s to 153.2 s (Table 2). These times enable to calculate the sterilization values corresponding to the holding phase:

𝐹 = 10 𝑇-𝑇 𝑟𝑒𝑓 𝑧 × 𝑡
(3) where F (unit of time) is the sterilization value, T (°C) is the treatment temperature, Tref (°C) is the reference temperature for sterilization (i.e. 121.1 °C), z (°C) is the temperature elevation that reduces the decimal reduction time by one log (i.e. 10 °C) and t (unit of time) is the duration of the treatment.

Experimental design

Two separate trials were done: first a resilience trial was dedicated to check the presence or absence of fouling in the ohmic cell. It consisted in running 1.5 h at a holding temperature of 120 °C and then increase temperature to 130 °C and 140 °C and run at 140 °C for 2.5 more hours. During the holding phase at 140 °C, samples were regularly taken at different residence times -S1, S4, S5, S7 (Table 2)and measured for FAST index and color for comparison with untreated infant formula.

The second trial was dedicated to the kinetic study, during which three holding temperatures were applied -120, 130 and 140 °Cand samples were taken at seven different heating times (Table 2) after stabilizing the system; each experiment was made in triplicate. All samples were analyzed for six chemical markers of the thermal treatment: FAST index, soluble proteins, color, vitamin C, furosine and CML.

Results and discussion

Process comparison between the two technologies and energetic considerations

The treatment of the product required enough energy to heat the infant formula from its initial temperature T0 (50 °C) until the target temperature Tf (120 -140 °C). This quantity of energy can be evaluated using Eq. ( 4) and knowing the flow rate, thermo-physical parameters (Table 1) and temperature difference between the initial and the final states of the product:

𝑄 ̇= 𝑚 IF ̇× 𝐶 p IF × (𝑇 f -𝑇 0 ) IF (4)
where 𝑄 ̇ (W) is the energy necessary to increase the temperature of the infant formula (IF) from its initial to its final temperature; 𝑚 IF ̇ (kg/s) is the mass flow of infant formula and 𝐶 p IF (J/kg/K) its specific heat capacity. Under such temperature conditions and with infant formula parameters calculated at a median temperature between T0 and Tf, Eq. ( 4) gives an energy value of 11 to 15 kW. This energy had to 9/16 be delivered either by electricity for the OH or by steam condensation and cooling for the SI. Knowing the temperature of the injected steam (Ts = 145 °C) it is then possible to calculate its flow rate for each final product temperature using Eq. ( 5):

𝑄 ̇= 𝑚 S ̇× ∆𝐻 (145 °C) + 𝑚 S ̇× 𝐶 p W × (𝑇 s -𝑇 f ) (5) where 𝑚 S ̇ (kg/s) is the mass flow of steam (or water); ∆𝐻 (180 °C) (J/kg) is the enthalpy of steam condensation at 145 °C; 𝐶 p W (J/kg/K) is the specific heat capacity of the water.

Eq. ( 5) gives steam flow rates of 0.005 and 0.007 kg/s for final product temperatures of 120 and 140 °C respectively, using the 𝑄 ̇ value given by Eq. ( 4) and ∆𝐻 (145 °C) and an estimated value of 𝐶 p W at a median temperature between Ts and Tf, according to [START_REF] Bergman | Fundamentals of heat and mass transfer[END_REF]. Such a steam flowrate allows predicting a non negligible dilution of the SI samples of 12 to 17%. The vaporization during sampling is expected to compensate, at least partly, the product dilution due to steam condensation. Furthermore, vaporization during sampling, although uncontrolled, is expected to be the same for OH and SI technology, which means that OH samples may be more concentrated than SI samples in this study, but to an unknown extent. Despite that, it was chosen not to correct the measured concentration of newly formed compounds (NFC) and to keep in mind that OH samples are probably slightly more concentrated than SI samples for accumulating NFC and slightly more diluted for transition markers. This uncertainty would be worth being removed in further studies.

The sterilization values (F) calculated for the holding phase show that the treatment becomes realistic compared to industrial standards (F > 152 s, corresponding to a 12 log reduction of a Clostridium botulinum population) for the samplings S2-S3 at 140 °C and S5 at 130 °C. For 120 °C, the treatment should have been continued a little longer after S7 (42-43 seconds).

Resilience trial

The resilience trial aimed at detecting any occurrence of fouling during the time covering the kinetic study. A first holding phase, set at 120 °C, lasted 1.5 hour and a second phase, set at 140 °C lasted 2.5 more hours. During the resilience trial, no temperature decay could be observed which indicates that no fouling occurred. Indeed, when a fouling layer appears on the electrodes, it acts like an insulation layer that provokes a loss in the quantity of energy delivered in the product and consequently temperature decay.

Besides, the use of chemical markers (FAST index and color) was also tested to check fouling during the duration of the 140 °C resilience test. Samples were taken from the S1, S4, S5 and S7 sampling points. The sampling point S1 (1.6 s) reflected the product state at the entrance of the holding section as well as the exit of the ohmic treatment cell (Fig. 1). The FAST index of S1 samples gave stable values around 20 during the 2.5 h running at 140 °C, above the value of 6.1 of the untreated sample. This indicates that no fouling occurred in the ohmic applicator while it most probably occurred in the holding section: the FAST index measured on S4 to S7 samples increased during the trial all along the holding tube (S4: 22.7 to 25.2; S7: 66.1 to 75.5). The color measurements were not sensitive enough to detect any change during the 2.5 h resilience trial at 140 °C. Such results confirm the good performances of this new ohmic applicator regarding fouling.

10/16

Determination of the kinetic parameters of the observed chemical markers

The evolution of the chemical markers concentration with time can be separated into two categories: the ones with linear increase of concentration during the isothermal phase of the thermal treatment (Fast index, color, furosine and CML) and those with constant concentration (vitamin C and soluble proteins). The first category regroups specific markers of Maillard reaction like furosine and CML while FAST index and color (divided in its 3 components L*, a* and b*) are global markers of the same reaction. The kinetic behavior of such markers has been successfully represented in a previous article by [START_REF] Roux | Kinetics of Maillard reactions in model infant formula during UHT treatment using a static batch ohmic heater[END_REF] using a model of pseudo-zero order (Eq. ( 7)), where k(T) parameter could be estimated using linear regression applied to the isothermal part of the curves (Table 3):

t k C C T    ) ( 0 (7)
C is the reactant concentration at time t and C0 its initial concentration; k(T) is the reaction rate depending on temperature T. Arrhenius parameters -kinetic rate constant k(Tref) at a reference temperature Tref chosen to be the middle of the studied temperature interval and activation energy Ea (Eq. ( 8))could then be estimated by fitting Eq. ( 8) on all experimental data. This reparametrized Arrhenius equation enables to be more accurate for the k(Tref) determination than for a determination of the preexponential factor (k0) classically used but very far away from the temperature domains generally studied [START_REF] Van Boekel | Kinetic modeling of reactions in foods[END_REF].

                     ref ref) ( ) ( 1 1 exp T T R E k k a T T (8)
Table 3. Kinetic parameters (mean ± standard deviation) of the different markers of Maillard reaction observed during the thermal treatment of an infant formula by ohmic heating or steam injection. k(T) were estimated using linear regression on the data obtained during the holding phase. k(Tref = 130 °C) and Ea were estimated by a global fitting of all experimental data according to Eq. ( 8).

Reaction rate of a pseudo-zero order type of kinetic (U/s) (a) Parameters calculated from Arrhenius equation k( 120 Comparable values of k(T), k(Tref) and Ea were obtained for steam injection and ohmic heating regarding FAST index. In the case of color (ΔL*, Δa* and Δb*), very little variations were detected all along the trials, resulting in very low kinetic rate constants that could not be well identified (high standard deviations). But the Arrhenius parameters showed no major differences between both technologies. Finally, furosine and CML presented quite high differences on the k values but activation energies were close. The results presented in Table 3 are fairly comparable to the kinetic data obtained in a previous work for a model infant formula treated by ohmic heating in a lab-scale batch reactor [START_REF] Roux | Kinetics of Maillard reactions in model infant formula during UHT treatment using a static batch ohmic heater[END_REF].

Comparison of the kinetic behavior of the markers between ohmic heating and steam injection

Thanks to the device presented in Fig. 1, the residence time of the liquid formula could be considered approximately the same for ohmic heating and steam injection. The thermal markers could thus be compared by plotting them for each technology, the closer to the bisector, the better the equivalence between both technologies (Fig. 3). Confirming the proximity of the kinetic parameters shown in Table 3, the results obtained on FAST index emphasized the equivalence of both technologies. No significant difference could be detected on soluble proteins content between ohmic heating and steam injection. Even though furosine and CML had comparable k0 and Ea, they seemed to be generated in higher amount by ohmic heating. However, this is only noticeable at high treatment times whereas industrial treatments would only last several seconds, i.e. in the equivalence zone. As regards color, slight differences could be detected in line with the kinetic parameters. Most of the L* points stayed on the zero position for ohmic heating while darkening was observed for steam injection. a* and b* (data not shown) were the same for both technologies except for extreme time/temperature conditions (sampling times 6 and 7 at 130 and 140 °C) where it showed higher values for ohmic heating than steam injection. Therefore, if we exclude the extreme holding times which are not realistic, almost no color differences were generated by both technologies and they would probably not be detectable on the final industrial product. Finally, vitamin C seemed to be better preserved by ohmic heating than steam injection: the data points are clearly situated above the bisector, showing significantly higher residual vitamin C after a continuous sterilization by ohmic heating than steam injection. However, such results should be confirmed because they can, at least partially, be due to the fact that OH samples are not diluted during the treatment contrary to SI ones.

12/16

Conclusion

Continuous ohmic heating applied at pilot scale gave very satisfactory results in terms of preservation of the nutritional quality of a liquid infant formula sterilized under UHT conditions. Chemical markers like FAST index, furosine or carboxymethyllysine indicated an equivalent extent of Maillard reaction during sterilization by ohmic heating in comparison with steam injection which is classically considered the best performing technology. The product sterilized by ohmic heating appeared to be slightly richer in Maillard products like furosine and CML than by steam injection but this only for extreme time/temperature settings which were not representative of real industrial conditions. The color seemed to remain more constant with ohmic heating than with steam injection. The extent of soluble protein denaturation was comparable for both technologies while vitamin C appeared to be even better preserved by ohmic heating. Keeping in mind that the markers concentrations for SI are probably slightly underestimated, the conclusion of equivalence between the two treatments makes even more sense. Furthermore, since no fouling was detected in the ohmic applicator after two and a half hours processing at 140 °C, continuous ohmic heating could thus be considered as a very promising technology for industrial thermal treatment of fragile dairy products. Such conclusions would of course deserve to be confirmed by trials with representative durations of treatment of the industrial conditions.

  1 ICARE : Impeding neo-formed Contaminants Accumulation to Reduce their health Effects; Collective research project(COLL-CT-2005-516415). Sixth Framework Program of the EuropeanCommunity (2002Community ( -2006) )

Fig. 1 .

 1 Fig. 1. Schematic diagram of the continuous sterilization device.

Fig. 2 .

 2 Fig. 2. Schematic diagram of the ohmic cell with coaxial symmetry according to the principle proposed in Pain et al. (2013).

Fig. 3 .

 3 Fig. 3. Comparison of FAST index, soluble proteins, furosine, carboxymethyllysine (CML), color (L*) and residual vitamin C of a liquid infant formula thermally treated under the same holding temperature conditions of sterilization by steam injection or ohmic heating; the full symbols correspond to predicted values of furosine and CML by front face fluorescence and the empty ones to chemical analysis by GC-MS-MS.

Table 1 .

 1 General composition and physical properties of the liquid infant formula.

	Composition	(g/L)
	Water	887
	Sugars	84
	Lactose		84
	Proteins	14
	Casein		11.5
	Whey		2.5
	Lipids	36
	Saturated fatty acids		14.6
	Mono-unsaturated fatty acids	12.8
	Poly-unsaturated fatty acids	7.1
	Minerals	4.6
	Iron		0.009
	Vitamins	0.43
	Vitamin C		0.00063
	Physical properties
	Density at 20 °C(a)

Table 2 .

 2 Average treatment times and sterilization values corresponding to the seven sampling points after ohmic heating or steam injection continuous sterilization of the liquid infant formula.

	Sampling point	Average treatment time (s) OH SI	120 °C	Sterilization values (F) (s) 130 °C OH SI	140 °C
	S1	1.6	0.4		12.4	3.1
	S2	2.0	0.8		15.5	6.2
	S3 S4 S5	5.9 13.7 39.5	4.8 12.5 38.3	𝐹 (130 °C) ÷ 10	45.8 106 307	37.3 97.0 297	𝐹 (130 °C) × 10
	S6	113.0	111.9		877	869
	S7	154.3	153.2		1198	1189

Acknowledgements

The authors would like to thank the European Community for their financial support to ICARE project (6 th Framework Programme, Grant No. COLL-CT-2005-516415) and the Stazione Sperimentale per l'Industria delle Conserve Alimentari (SSICA) in Parma Italy, for housing the experimental part of this work.