Alexandre Ern

Jean-Luc Guermond

Weighting the edge stabilization

Keywords: Finite elements, conservation equations, edge stabilization, linear stabilization, nonlinear viscosity. AMS subject classifications. 65M12, 35L65, 65M60

A modification of the edge stabilization technique is proposed to improve the behavior of the method when solving conservation equations with non-smooth data and/or non-smooth solutions. The key ingredient is tempering the edge stabilization in regions of large gradients through appropriate weights. The new method is shown to preserve the convergence properties of the original method on smooth solutions and numerical tests indicate that it performs better on non-smooth solutions.

1. Introduction. Linear stabilization techniques are known to be very effective for solving linear first-order PDEs like the transport equation. In particular, denoting h the meshsize and k the polynomial degree of the approximation, linear stabilization techniques yield the near optimal convergence rate O(h k+ 1 2 ) in the L 2 -norm for smooth solutions. The situation is not so bright when it comes to solving linear problems with non-smooth data and nonlinear conservations equations with non-unique weak solutions. Linear stabilization methods generally promote the Gibbs phenomenon and introduce high-order dissipation that, in the case of nonlinear conservation equations, can lead to convergence to non-entropic solutions.

We focus our attention on the edge stabilization technique [START_REF] Burman | A unified analysis for conforming and nonconforming stabilized finite element methods using interior penalty[END_REF][START_REF] Burman | Edge stabilization for Galerkin approximations of convection-diffusion-reaction problems[END_REF] as a prototype of linear stabilization, that is relatively easy to implement with H 1 -conforming finite elements. We show through numerical examples that edge stabilization promotes the Gibbs phenomenon, and, for nonlinear conservation equations with non-convex flux, cannot select the proper entropic solution. We also show that, when not properly scaled, the extra dissipation induced by edge stabilization can transform a convergent method into a non-convergent one. The purpose of this paper is then to introduce and analyze a modified version of edge stabilization that does not suffer from the above problems. The main modification is tempering, through appropriate weights, the edge stabilization in regions where the discrete solution exhibits large gradients. This may seem a bit counter-intuitive at first glance, since the use of linear stabilization techniques is often motivated to counter spurious oscillations that are produced by large gradients. The proposed method is proved to deliver the near optimal convergence rate O(h k+ 1 2 ) in the L 2 -norm for smooth solutions. Numerical tests on the linear transport equation in one and two space dimensions show that the weighted edge stabilization performs as required when combined with a nonlinear viscosity method: it no longer promotes the Gibbs phenomenon and does not prevent the nonlinear viscosity method to converge to the correct entropic solution. In other words, the weighted edge stabilization does not antagonize the nonlinear viscosity method. Quite importantly, when combined with a nonlinear viscosity method and for polynomial orders larger than or equal to two, we observe that the weighted edge stabilization increases the convergence order of the nonlinear method in the regions where the solution is smooth. Thus, when combining the weighted edge stabilization with a nonlinear viscosity method, one improves the convergence order of the nonlinear viscosity method without sacrificing its weakened maximum principle property and its ability to properly converge to entropic solutions.

The paper is organized as follows. In §2, we set the notation and present numerical experiments illustrating the main difficulties that are addressed herein. In §3, we introduce and analyze the weighted edge stabilization method. In §4, we present one-and two-dimensional tests to illustrate the improvements achieved by weighting the edge stabilization. Finally, we draw some conclusions in §5.

2.

Preliminaries. The objective of this section is twofold: (i) to set the notation and the model problems we are interested in; (ii) to present numerical experiments that identify the main difficulties that we want address in the present work.

2.1. Formulation of the problem. We are interested in approximating the solution of scalarvalued conservation equations in the form (2.1)

∂ t u + ∇•f (u) = 0, u(x, 0) = u 0 (x), (x, t) ∈ Ω×R + where Ω is an open polyhedral domain in R d and f ∈ C 1 (R; R d ).
For the sake of simplicity, we assume that there are no issues with the boundary conditions; for instance, either we assume periodic boundary conditions or the initial data is compactly supported and we are interested in the solution before the domain of dependence of u 0 reaches the boundary of Ω. We assume that (2.1) has a unique entropic solution satisfying additional entropy inequalities that we omit for brevity.

In order to approximate the entropic solution of (2.1) with H 1 -conforming finite elements, we consider a mesh family {K h } h>0 that we assume to be conforming (no hanging nodes) and shaperegular in the sense of Ciarlet. By convention, the elements in {K h } h>0 are closed in R d . The reference element is denoted K and the map between K and an arbitrary element K ∈ K h is denoted Φ K : K -→ K. We define the scalar-valued finite element approximation space (2.2)

X h = {v ∈ C 0 (Ω; R); v| K •Φ K ∈ P k , ∀K ∈ K h },
where k ∈ N\{0} and P k denotes the set of multivariate polynomials of total degree at most k. For all K ∈ K h , h K denotes the diameter of K divided by k and h := max K∈K h h K is the so-called meshsize.

Let t F > 0 be the final simulation time. We call Galerkin solution of (2.1) the function

u h ∈ C 1 ([0, t F ]; X h ) such that (2.3) Ω v∂ t u h dΩ + Ω v∇•f (u h ) dΩ = 0,
for all v ∈ X h and all t ∈ (0, t F ), and u h (t = 0) = u 0,h , where u 0,h is an appropriate approximation of u 0 in X h . In all the numerical tests reported herein, the time stepping is done either with the so-called SSP RK3 method, see e.g., [START_REF] Gottlieb | Strong stability-preserving high-order time discretization methods[END_REF], or the standard RK4 method. To avoid mixing the space and time discretization errors, we always perform our tests with a small CFL, say CFL = 0.2 or even less. We finally emphasize that the mass matrix is never lumped in our simulations.

The Galerkin solution is known to be a poor approximation of the solution to (2.1) even when the flux is linear. One simple device to stabilize the approximation is to add first-order viscous dissipation. Henceforth, we call viscous solution the function

u h ∈ C 1 ([0, t F ]; X h ) such that u h (t = 0) = u 0,h and (2.4) Ω v∂ t u h dΩ + Ω v∇•f (u h ) dΩ + n visc (u h ; v) = 0, for all v ∈ X h and all t ∈ (0, t F ), with (2.5) n visc (w; v) := c max K∈K h h K f ′ (w) L ∞ (K) K ∇w•∇v dK.
We take c max = 1 2k in one space dimension and c max = 1 4k in two space dimensions on triangular meshes.

It is well-known that the viscous solution is only first-order accurate. The performance of the method can be greatly improved by substituting the first-order viscous dissipation by some linear and/or nonlinear stabilization mechanism. We are going to use in this paper the so-called entropy viscosity [START_REF] Guermond | Entropy-based nonlinear viscosity for Fourier approximations of conservation laws[END_REF][START_REF] Guermond | On the use of the notion of suitable weak solutions in CFD[END_REF][START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF] as a nonlinear stabilization mechanism. We call entropy viscosity solution the function

u h ∈ C 1 ([0, t F ]; X h ) such that u h (t = 0) = u 0,h and (2.6) Ω v∂ t u h dΩ + Ω v∇•f (u h ) dΩ + n ev (u h ; v) = 0, for all v ∈ X h and all t ∈ (0, t F ), with (2.7) n ev (w; v) := K∈K h ν K (w) K ∇w•∇v dK,
where ν K (•) : C 1 ([0, t F ]; X h ) -→ P 0 (K) is a nonlinear viscosity functional. The details on how the nonlinear viscosity is constructed are reported in the Appendix. At this point it is not important to go through the detailed construction of n ev (•; •); it suffices to know that the entropy viscosity solution has reasonable convergence properties. For instance, the entropy viscosity solution has been observed to satisfy a weakened maximum principle in the sense that, for all ǫ > 0, there is h 0 > 0 and there are uniform constants c and α > 0 such that

(2.8) u h (t) L ∞ (Ω) ≤ u 0 L ∞ (Ω) + ch α , ∀h < h 0 , ∀t > ǫ.
This property is illustrated, for instance, in Table 2.1(a) below. When f (u) = βu, with R d -valued velocity field β, (2.1) reduces to the linear transport equation. Many linear (symmetric) stabilization techniques can be considered for solving the linear transport equation, e.g.,, subgrid viscosity [START_REF] Guermond | Stabilization of Galerkin approximations of transport equations by subgrid modeling[END_REF], edge stabilization [START_REF] Burman | A unified analysis for conforming and nonconforming stabilized finite element methods using interior penalty[END_REF][START_REF] Burman | Edge stabilization for Galerkin approximations of convection-diffusion-reaction problems[END_REF][START_REF] Burman | Explicit Runge-Kutta schemes and finite elements with symmetric stabilization for first-order linear PDE systems[END_REF], and discontinuous Galerkin [START_REF] Johnson | An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF][START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF][START_REF] Hesthaven | Nodal discontinuous Galerkin methods[END_REF][START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]. All these methods are known to yield the near optimal convergence rate O(h k+ 1 2 ) in the L 2 -norm for smooth solutions. We focus herein on the edge stabilization technique, which we think is relatively simple to implement with H 1 -conforming finite elements. Let F i h be the set of all mesh interfaces. By convention, interfaces are closed in R d-1 . For all F ∈ F i h , we denote h F the diameter of F . Denoting K 1,F ∈ K h , K 2,F ∈ K h the two (closed) cells so that F = K 1,F ∩ K 2,F , we set ∆ F = K 1,F ∪ K 2,F , see Figure 3.1 below. Moreover, letting v be a scalar-valued function defined over ∆ F and continuous over K 1,F and K 2,F , we denote {v} (x) = 1 2 (v| K1,F (x) + v| K2,F (x)) the average of v at F . The edge stabilization bilinear form is defined as follows:

(2.9)

n ed (w; v) = c ed F ∈F i h h 2 F f ′ (w) L ∞ (F ) F {∂ n w} {∂ n v} dF,
where ∂ n is the outward normal derivative (so that {∂ n w} is proportional to the jump of the normal gradient across F ) and c ed is a user-defined parameter. In our numerical experiments, we set c ed = 0.05. We call edge stabilized solution the function

u h ∈ C 1 ([0, t F ]; X h ) such that u h (t = 0) = u 0,h and
(2.10)

Ω v∂ t u h dΩ + Ω v∇•f (u h ) dΩ + n ed (u h ; v) = 0.
for all v ∈ X h and all t ∈ (0, t F ). The objective of this paper is to investigate some aspects of the stabilization properties of the bilinear form n ed (•; •). We want to show that, although this stabilization technique performs extremely well in smooth regions, it has counter-productive effects in regions of shocks and large gradients. The purpose of the rest of this section is to illustrate through numerical experiments some of the negative effects of n ed (•; •). A weighting technique that cures these problems is proposed and analyzed in §3.

2.2. One-dimensional transport. We consider the one-dimensional transport problem (2.11)

∂ t u + ∂ x u = 0, u(x, 0) = 1, if 0.4 < x < 0.7, 0, otherwise,
over the interval Ω = (0, 1) with periodic boundary conditions. We compute the solution at t F = 1 with continuous P 1 finite elements on various uniform meshes. To assess departures from the maximum principle, we compute the following indicators:

(2.12)

e Max := max x∈Ω u h (x, 1) -1, e Min := -min x∈Ω u h (x, 1
).

Edge stabilization alone.

We first compute the Galerkin solution on five uniform meshes composed of 100, 2×100, . . . , 2 4 ×100 cells and we set CFL = 0.2. The graphs of the solutions are shown in Figure 2.1(a). We observe the familiar spurious oscillations that characterize the Galerkin technique. The graphs of the solutions computed with edge stabilization are shown in Figure 2.1(b). This test exemplifies at the same time the stabilizing capability of the edge stabilization and its inability to counter the so-called Gibbs phenomenon triggered in the vicinity of discontinuities and large gradients. Figure 2.1(c) displays details of the graphs of the solutions in the region x ∈ [0.65, 0.71]. Numerical tests on eight refinement levels (nor reported here) show that both indicators e Max and e Min are bounded away from 0, i.e., mesh refinement does not to help satisfy the maximum principle.

Edge stabilization plus entropy viscosity.

It is frequently advocated in the literature that linear stabilization must be supplemented with a shock capturing technique to handle properly shocks and large gradients [START_REF] Hughes | A new finite element formulation for computational fluid dynamics. IV: A discontinuity-capturing operator for multidimensional advective-diffusive systems[END_REF][START_REF] Johnson | On the convergence of a finite element method for a nonlinear hyperbolic conservation law[END_REF][START_REF] Codina | A discontinuity-capturing crosswind-dissipation for the finite element solution of the convectiondiffusion equation[END_REF][START_REF] Burman | Nonlinear diffusion and discrete maximum principle for stabilized Galerkin approximations of the convection-diffusion-reaction equation[END_REF][START_REF] Guermond | Subgrid stabilized projection method for 2D unsteady flows at high Reynolds number[END_REF]. We want to investigate the effects of combining the edge stabilization with the entropy viscosity described in the Appendix.

We perform the following tests. We compute the entropy viscosity solution and the edgestabilized entropy viscosity solution of the one-dimensional transport equation (2.11) using P 1 finite The results for the entropy viscosity solution and the edge-stabilized entropy viscosity solution are reported in Table 2.1(a) and Table 2.1(b), respectively. We observe that both indicators e Max and e Min for the entropy viscosity solution converge to zero with the meshsize, in agreement with the claim made in Section 2.1 that the entropy viscosity solution satisfies a weakened maximum principle in the form (2.8). On the other hand, we observe in Table 2.1(b) that the weakened maximum principle is lost when edge stabilization is added to the entropy viscosity. Thus, by adding edge stabilization to a method that satisfies a weakened maximum principle, we obtain a method that does not satisfy the maximum principle, even in the weak sense defined above.

The above problem can be easily fixed by increasing the strength of the entropy viscosity. For instance, the weakened maximum principle can be recovered by using γ = 1 2 in the following definition of the entropy viscosity

(2.13) ν K (w) = min(c max h K f ′ (w) L ∞ (K) , c ev h γ K R K (w))
, where R K (•) is the entropy residual (see the Appendix). The definition of ν K introduced in the Appendix uses γ = 1, which makes it asymptotically smaller by a factor h 1 2 than the above definition. The γ = 1 2 fix is marginally satisfactory since it makes the method more dissipative and deteriorates its convergence properties. For instance, convergence tests on the one-dimensional transport problem (2.11) with P 1 finite elements reveal that the convergence rate of the entropy viscosity method in the L 1 -norm is 3 4 with γ = 1 and 2 3 with γ = 1 2 . 2.3. One-dimensional non-convex conservation equation. We now consider a problem with non-convex flux proposed in [START_REF] Kurganov | Adaptive semidiscrete central-upwind schemes for nonconvex hyperbolic conservation laws[END_REF] to test the edge stabilization technique with nonlinear conservation equations. We restrict ourselves to the one-dimensional domain Ω = (0, 1) and we consider the following scalar flux and initial data:

(2.14) f (u) = 1 4 u(1 -u) if u < 1 2 , 1 2 u(u -1) + 3 16 if 1 2 ≤ u, u 0 (x) = 0, x ∈ [0, 0.35], 1, x ∈ (0.35, 1].
The entropic solution to this problem is a composite wave composed of a shock followed by a rarefaction wave. This problem is challenging since many second-order central schemes with compressive limiters are known to converge to weak solutions that are not entropic. For instance, it is demonstrated in [START_REF] Kurganov | Adaptive semidiscrete central-upwind schemes for nonconvex hyperbolic conservation laws[END_REF] that the so-called central-upwind scheme using second-order piecewise linear reconstruction with either the superbee limiter or the so-called minmod2 limiter fails to converge to the entropic solution. We compute the solution at time t F = 1 with continuous P 1 finite elements. The solution at t F = 1 is composed of a shock wave located at x s (1) = 1 4 ( √ 6 -1) followed by a rarefaction wave. The left limit of the solution at the shock is u - s (t F ) = 0 and the right limit

is u + s (t F ) = 2 x1(tF)-xs(tF) (x 1 (t F ) -x 0 -3 16 t F ) -1
, where x 0 = 0.35 and x 1 (t) = 1 2 t + x 0 is the time-dependent location of the head of the rarefaction wave.

2.3.1. Edge stabilization alone. We show in Figure 2.2 the Galerkin and the edge stabilized solutions obtained at t F = 1 on a uniform mesh composed of 1000 cells. The time stepping is done with the SSP RK3 scheme, and to avoid time discretization errors, the time step size is based on CFL = 0.01. The solution shown in the left panel of Figure 2.2 is the Galerkin solution and that shown in the right panel is the edge stabilized solution. The entropic solution is shown in blue solid line. It is clear that none of these approximations converge to the entropic solution in any possible norm. The edge stabilized solution is almost free of spurious oscillations but converges to a non-entropic weak solution. 2.3.2. Edge stabilization plus entropy viscosity. We show in Figure 2.3 the entropy viscosity solution and the entropy viscosity solution with edge stabilization at t F = 1 on a uniform mesh composed of 1000 cells. The entropy viscosity solution is shown in the left panel and the entropy viscosity solution with edge stabilization is shown in the center and right panels. The exact solution (so-called entropic solution) is shown in blue solid line. This tests show that the entropy viscosity solution converges to the entropic solution, whereas the entropy viscosity solution with edge stabilization does not. Thus, by adding edge stabilization to a method that converges to the correct weak solution, we obtain a method that converges to a wrong weak solution. This result is similar to what has been observed in [START_REF] Kurganov | Adaptive semidiscrete central-upwind schemes for nonconvex hyperbolic conservation laws[END_REF] concerning the second-order piecewise linear reconstruction combined with either the superbee or the minmod2 limiter.

Edge stabilization plus first-order viscosity.

In this section, we present two numerical examples showing that the edge stabilization can have adverse effects even on the first-order viscosity method. We first consider the one-dimensional inviscid Burgers equation

(2.15) ∂ t u + ∂ x ( 1 2 u 2 ) = 0, u(x, 0) = sin(2πx),
over the interval Ω = (0, 1) with periodic boundary conditions. The solution is computed at t F = 0.25 with continuous P 1 finite elements on a mesh composed of 200 elements, and CFL = 0.025. We compute the first-order viscous solution and the first-order viscous solution with edge stabilization. The coefficient c ed in (2.9) is set to 1, and the coefficient c max in (2.5) to 1 2 . When using finite differences on a one-dimensional uniform grid, setting c max = 1 2 corresponds to replacing the centered differences by first-order upwind differences, and the resulting scheme is known to be monotone. The graph of the two solutions is shown in Figure 2.4(a). For the solution with edge stabilization, we observe over-shoots and under-shoots, and the amplitude of these spurious features is constant as the mesh is refined. This is characteristic of the Gibbs phenomenon. These over-shoots and undershoots can be tamed by increasing the viscous dissipation beyond what should normally be necessary. Numerical tests (not shown) reveal that c max = 2 is the lower bound on c max that makes the viscous dissipation strong enough to overcome the Gibbs phenomenon. These results indicate that edge stabilization (and possibly linear stabilization at large) tends to promote the Gibbs phenomenon. The adverse effects of edge stabilization are even more dramatic on conservation equations with non-convex flux. We consider again the test case of §2.3. We show in Figure 2.4(b) a zoom of the graph of the solutions obtained with first-order viscosity (c max = 1 2 ) without edge stabilization and with edge stabilization (c ed = 1) computed on two uniform meshes composed of 4, 000 and 10, 000 cells. We observe that the viscous solution converges to the entropic solution (as expected), whereas the edge stabilized solution converges to a plateau between the expansion wave and the shock, which is wrong. The same effect is observed in Figure 2.4(c) with the choice c ed = 0.05 in (2.16)

n ed (w; v) = c ed f ′ (w) L ∞ (Ω) F ∈F i h h 2 F F
{∂ n w} {∂ n v} dF, 2.5. Conclusions from numerical tests. The first series of tests presented in §2.2.1 shows that the linear edge stabilization does a great job at suppressing spurious oscillations in the regions where the solution is smooth, but this technique cannot get rid of the Gibbs phenomenon. The second series of tests reported in §2.2.2 show that the linear edge stabilization actually promotes the Gibbs phenomenon. The third series of tests on the one-dimensional nonlinear scalar conservation equation with a non-convex flux in §2.3.1 demonstrates that the linear edge stabilization does again a great job at removing the spurious oscillations plaguing the Galerkin solution, but it does not have the correct type of dissipation to make the approximate solution converge to the entropic solution. Finally, the tests reported in §2.3.2 and §2.4 show that not only the linear edge stabilization does not produce the right dissipation, but the type of dissipation that it produces can transform a convergent method (either first-order linear viscosity or nonlinear entropy viscosity) into a non-convergent one. The authors conjecture that the above conclusions are not restricted to edge stabilization but can be extended to some (or most of the) other linear stabilization methods available in the literature.

3.

Weighting the edge stabilization. We introduce in this section a weighting technique for the edge stabilization, and we prove that its convergence properties on the linear transport problem are identical to that of the original unweighted edge stabilization in the case of smooth solutions. We consider the linear transport equation

(3.1) ∂ t u + ∇•(βu) = 0, u(x, 0) = u 0 (x), (x, t) ∈ Ω×R + ,
in space dimension d = 2 or d = 3. We assume that β is Lipschitz and divergence-free in Ω. Recall that for simplicity, we assume either periodic boundary conditions, compactly supported solutions, or

β•n| ∂Ω = 0.
In what follows, a b means that inequality a ≤ c b holds with a constant c independent of h (but possibly depending on the mesh-regularity, the polynomial degree k, and the regularity of the problem data and the exact solution). Without loss of generality, we assume h ≤ 1. For any set R ⊂ Ω (a mesh element, a mesh face, or a collection thereof), we denote by • L p (R) the usual L p (R)-norm, 1 ≤ p ≤ ∞, for scalar-or vector-valued functions.

Principle of the method. Consider the following discrete solution u

h ∈ C 1 ([0, t F ]; X h ), with X h defined in (2.2), so that u h (t = 0) = u 0,h and (3.2) Ω v∂ t u h dΩ + Ω v∇•(βu h ) dΩ + n lim,ed (u h ; u h , v) = 0.
for all v ∈ X h and all t ∈ (0, t F ). The weighted edge stabilization semi-linear form is defined by

(3.3) n lim,ed (z; w, v) = c ed F ∈F i h α(g F (z))h 2 F |β| F F {∂ n w} {∂ n v} dF,
with the shorthand notation |β| F := β L ∞ (F ) and g F (z) is a measure of the gradient of z around F which we take in the form

(3.4) g F (z) = ℓ -1 | ∇z ∆F |,
where φ R := meas(R) -1 R φ dR denotes the average of a function φ over a set R ⊂ Ω and where the global scaling parameter ℓ is, e.g., set to ℓ := | ∇u 0,h Ω |. The key ingredient in (3.3) is the function α : R + → (0, 1] which weights the amount of edge stabilization. The function α must be such that (i) α is non-increasing; (ii) there is α 0 ∈ R + and λ ∈ R + such that, for all r ≥ 1, α(r) ≥ α 0 r -λ . Typically, α(r) → 0 as r → +∞ so as to turn off edge stabilization in regions of large gradients. Condition (ii) then means that α must not decrease too quickly, so as to retain the optimal convergence properties of the method for smooth solutions.

Convergence analysis.

We first prove a convergence result in two dimensions. In this situation, there is no restriction on the parameter λ controlling the decrease of the weighting function α at infinity. Theorem 3.1. Let u and u h be the solutions to (3.1) and to (3.2), respectively. Assume that

u ∈ C 1 (0, t F ; H k+1 (Ω)) ∩ C 0 (0, t F ; W k+1,∞ (Ω)). Assume d = 2 and quasi-uniform meshes. Let λ > 0. Then, for all t ∈ [0, t F ], (3.5) (u -u h )(•, t) L 2 (Ω) + t 0 n lim,ed (u h ; u h , u h ) dτ 1 2 h k+ 1 2 .
Proof. The proof is decomposed into three steps.

Step 1: error equation. For all t ∈ [0, t F ], let w(•, t) be the L 2 -orthogonal projection of the exact solution u(•, t) onto the discrete space X h (in the case of compactly supported solutions, we project onto X h ∩ H 1 0 (Ω)). We define the quantities (the dependence with respect to t is now left implicit)

e := u h -w, η := u -w, so that the approximation error is u h -u = e -η. Observing that Ω v∂ t w dΩ + Ω v∇•(βw) dΩ = - Ω v∂ t η dΩ - Ω v∇•(βη) dΩ, ∀v ∈ X h ,
=: T 1 (η, v) + T 2 (η, v) + T 3 (u h ; η, v), (3.6)
where we have used the fact that, for the exact solution u, {∂ n u} = 0 for all F ∈ F i h , so that

n lim,ed (u h ; u h , v) = n lim,ed (u h ; e, v) -n lim,ed (u h ; η, v) + n lim,ed (u h ; u, v)
= n lim,ed (u h ; e, v)n lim,ed (u h ; η, v).

Step 2: basic estimates. Testing (3.6) with v = e and using the conservativity property

Ω e∇•(βe) dΩ = 0,
which holds owing to the choice of boundary conditions and the fact that β is divergence-free, we infer that

1 2 d dt e 2 L 2 (Ω) + n lim,ed (u h ; e, e) ≤ |T 1 (η, e)| + |T 2 (η, e)| + |T 3 (u h ; η, e)|.
Moreover, since u ∈ C 1 (0, t F ; H k+1 (Ω)), classical finite element interpolation properties [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF] yield

(3.7) h K ∇η L 2 (K) + η L 2 (K) + ∂ t η L 2 (K) h k+1 K , ∀K ∈ K h .
Hence, by the Cauchy-Schwarz inequality, we obtain

|T 1 (η, e)| h k+1 e L 2 (Ω) , |T 2 (η, e)| h k e L 2 (Ω) , |T 3 (u h ; η, e)| ≤ n lim,ed (u h ; e, e) 1 2 n lim,ed (u h ; η, η) 1 2 h k+ 1 2 n lim,ed (u h ; e, e) 1 2
,

where the last bound results from the fact that n lim,ed (u h ; η, η)

1 2 ≤ n ed (η, η) 1 2
(since α F (u h ) ≤ 1) and the classical bound n ed (η, η)

1 2 h k+ 1 2 .
Collecting the above estimates, using Young's inequality, and recalling that h ≤ 1, we arrive at

d dt e 2 L 2 (Ω) + n lim,ed (u h ; e, e) h k e L 2 (Ω) + h 2k+1 ,
whence by using Gronwall's Lemma together with h ≤ 1, we obtain the suboptimal estimate

(3.8) ∀t ∈ [0, t F ], e(•, t) L 2 (Ω) + t 0 n lim,ed (u h ; e, e) dτ 1 2 h k .
Step 3: improved estimate on T 2 (η, e). Let ǫ ≥ 0 and c 0 ≥ 1 (the value of these quantities is chosen later on). We fix a time t ∈ (0, t F ). We define the sets

F ♯ h := {F ∈ F i h ; g F (u h ) ≥ c 0 h -ǫ }, K ♭ h := {K ∈ K h ; ∀F ∈ F K , F ∈ F ♯ h }, K ♯ h := {K ∈ K h ; ∃F ∈ F K , F ∈ F ♯ h } = K h \ K ♭ h ,
where F K is the collection of all the mesh interfaces having a nonempty intersection with K (see Figure 3.2). Let β h be the continuous, piecewise affine interpolant of β on the mesh K h and set z h := β h •∇e.

F

Observe that z h is a piecewise polynomial of degree ≤ k, but it does not belong to X h because it can jump across interfaces. We define I av (z h ) ∈ X h to be so that its value at any Lagrange node is equal to the arithmetic average of the values taken by z h from the mesh elements sharing that node. It is well-known, see [START_REF] Achdou | A priori and a posteriori analysis of finite volume discretizations of Darcy's equations[END_REF][START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems[END_REF][START_REF] Burman | Continuous interior penalty hp-finite element methods for advection and advectiondiffusion equations[END_REF], that this type of averaging leads to the estimate, for all K ∈ K h , (3.9)

z h -I av (z h ) L 2 (K) F ∈FK h 1 2 F [[z h ]] L 2 (F ) ≤ F ∈FK h 1 2 F |β| F {∂ n e} L 2 (F ) ,
where [[•]] denotes the jump across F and where we have used the continuity of β h across F . Integrating by parts, accounting for boundary conditions, and using the fact that η is L 2 -orthogonal to X h , T 2 (η, e) can be decomposed as

T 2 (η, e) = - Ω η(β•∇e) dΩ = - Ω η(β -β h )•∇e dΩ - Ω ηz h dΩ = - Ω η(β -β h )•∇e dΩ - Ω ♭ η(z h -I av (z h )) dΩ - Ω ♯ η(z h -I av (z h )) dΩ =: T 2,0 (η, e) + T 2,1 (η, e) + T 2,2 (η, e),
where

Ω ♭ := ∪ K∈K ♭ h K and Ω ♯ := ∪ K∈K ♯ h K.
Since β is Lipschitz, the term T 2,0 (η, e) is readily bounded as follows:

|T 2,0 (η, e)| η L 2 (Ω) h ∇e L 2 (Ω) η L 2 (Ω) e L 2 (Ω) h 2k+1 ,
where we have used an inverse inequality and the estimate (3.8) on e L 2 (Ω) . For the term T 2,1 (η, e), we use (3.9), the Cauchy-Schwarz inequality, h K h F , and |β| F 1 to infer

|T 2,1 (η, e)|   K∈K ♭ h [α(g K (u h ))h K ] -1 η 2 L 2 (K)   1 2   K∈K ♭ h α(g K (u h )) F ∈FK h 2 F |β| F {∂ n e} 2 L 2 (F )   1 2 , with g K (u h ) := max F ∈FK g F (u h ). Since the function α is non-increasing, α(g K (u h )) ≤ α(g F (u h ))
for all F ∈ F K , so that this estimate implies

|T 2,1 (η, e)|   K∈K ♭ h [α(g K (u h ))h K ] -1 η 2 L 2 (K)   1 2 n lim,ed (u h ; e, e) 1 2
.

We now use Assumption (ii) on the function α to infer that, for all K ∈ K ♭ h and all F ∈ F K , α(g F (u h )) ≥ α(c 0 h -ǫ ) ≥ α 0 c -λ 0 h ǫλ since c 0 h -ǫ ≥ 1 (because c 0 ≥ 1, ǫ ≥ 0, and h ≤ 1). Thus, we obtain α(g

K (u h )) ≥ α 0 c -λ 0 h ǫλ , whence   K∈K ♭ h [α(g K (u h ))h K ] -1 η 2 L 2 (K)   1 2 h k+ 1 2 -1 2 λǫ , so that |T 2,1 (η, e)| h k+ 1 2 -1 2 λǫ n lim,ed (u h ; e, e) 1 2
.

We now bound T 2,2 (η, e). Since η L ∞ (Ω) h k+1 owing to u ∈ W k+1,∞ (Ω) and the approximation properties in L ∞ (Ω) of the L 2 -orthogonal projection for d = 2 and quasi-uniform meshes, see [START_REF] Crouzeix | The stability in Lp and W 1 p of the L 2 -projection onto finite element function spaces[END_REF], we infer

|T 2,2 (η, e)| ≤ η L ∞ (Ω ♯ ) meas(Ω ♯ ) 1 2 z h -I av (z h ) L 2 (Ω ♯ ) h 2k meas(Ω ♯ ) 1 2 ,
where we have used

z h -I av (z h ) L 2 (Ω) z h L 2 (Ω) ∇e L 2 (Ω) h k-1
. We now estimate meas(Ω ♯ ). This is achieved by proving meas(Ω ♯ ) meas(Ω ∆♯ ) and estimating meas(Ω ∆♯ ). Here, Ω ∆♯ := ∪ F ∈F ♯ h ∆ F , so that Ω ∆♯ is the union of all the elements that have a face in F ♯ h . Note that Ω ∆♯ ⊂ Ω ♯ ; see Figure 3 Consider K ∈ K ♯ h . By definition, there is F ∈ F K such that F ∈ F ♯ h . By definition also, ∆ F ⊂ Ω ∆♯ and K ∩ ∆ F = ∅ (recall that K and ∆ F are closed). This means that there is

K ′ ∈ Ω ∆♯ such that K ∩ K ′ is nonempty, so that K ♯ h ⊂ ∪ K ′ ∈Ω ∆♯ {K ′′ ∈ K h ; K ′ ∩ K ′′ = ∅}.
Hence, using mesh regularity, we infer meas(Ω ♯ ) ≤

K ′ ∈Ω ∆♯ meas({K ′′ ∈ K h ; K ′ ∩ K ′′ = ∅}) K ′ ∈Ω ∆♯ meas(K ′ ) = meas(Ω ∆♯ ).

The definition of F

♯ h now implies that meas(Ω ∆♯ )h -2ǫ ≤ F ∈F ♯ h meas(∆ F )h -2ǫ ≤ (c 0 ℓ) -2 F ∈F ♯ h meas(∆ F )| ∇u h ∆F | 2 ≤ (c 0 ℓ) -2 F ∈F ♯ h ∆F |∇u h | 2 dΩ ≤ 2(c 0 ℓ) -2 ∇u h 2 L 2 (Ω ∆♯ ) ≤ 6(c 0 ℓ) -2 ( ∇e 2 L 2 (Ω ∆♯ ) + ∇η 2 L 2 (Ω ∆♯ ) + ∇u 2 L 2 (Ω ∆♯ ) ).
Since ∇e L 2 (Ω ∆♯ ) h k-1 , ∇η L 2 (Ω ∆♯ ) h k , and ∇u L 2 (Ω ∆♯ ) ≤ meas(Ω ∆♯ ) 1 2 ∇u L ∞ (Ω ∆♯ ) , choosing c 0 = max( √ 12ℓ -1 ∇u L ∞ (Ω×(0,tF)) , 1) allows one to hide ∇u 2 L 2 (Ω ∆♯ ) on the left-hand side, so that meas(Ω ∆♯ ) h 2(k-1+ǫ) , and hence

|T 2,2 (η, e)| h 3k-1+ǫ .
Collecting the above estimates and dropping higher-order terms, we arrive at d dt e 2 L 2 (Ω) + n lim,ed (u h ; e, e) h 3k-1+ǫ + h 2k+1-λǫ .

For k ≥ 2, we can take ǫ = 0 and obtain the desired estimate of order h 2k+1 since 3k -1 ≥ 2k + 1 in this case. In the case k = 1, the two terms on the right-hand side are balanced by choosing ǫ = 1 1+λ . This leads to the estimate (3.8) with the sharper bound h 1+ρ with ρ = 1 2 ǫ. We now use a bootstrap argument. Suppose that estimate (3.8) holds with bound h 1+ρn . Then, proceeding as above with a parameter ǫ n to be chosen, the new bound for T 2,2 (η, e) becomes h 2+ǫn+2ρn , while the bound for T 2,1 (η, e) is still h 3-λǫn . Balancing the two bounds leads to the choice

ǫ n = 1 1+λ (1 -2ρ n ), thus improving estimate (3.8) to h 1+ρn+1 with ρ n+1 = 1 2 (1 -λǫ n ) = 1 2 ( 1 1+λ + 2λ 1+λ ρ n ).
This recursive relation shows that ρ n converges to 1 2 (for all λ), thereby completing the proof. and initial data specified by (2.14). The computations are performed on five uniform meshes composed of 100, 200, 400, 800, and 1600 cells using the entropy viscosity method with weighted edge stabilization. The graphs of the solutions obtained at t F = 1 are shown in Figure 4.1. We observe that, contrary to what is shown in Figure 2.2(b), the method now converges to the correct entropic solution thereby confirming again that weighting the edge stabilization is a good idea. layer between the shock and the expansion wave. We have verified by mesh refinement that this approximation does not converge to the entropic solution. Finally, the graph of the solution with weighted edge stabilization is shown in Figure 4.2(c). We observe that the unphysical layer has been pushed back to the shock, thereby recovering the correct physical behavior.

Conclusions.

To conclude we would like to offer some thoughts on the relative importance of nonlinear viscosities (also called shock capturing in the literature) and linear stabilization techniques. Substantial efforts have been devoted to the construction of linear stabilization techniques (GaLS, RFB, SUPG, VMS, Streamline Diffusion, Subgrid Viscosity, Edge Stabilization, local projection, etc.), and it is often implicitly understood that linear stabilization is the workhorse whereas shock capturing is only meant to remove remaining spurious oscillations. As a result, the amount of research dedicated to constructing and analyzing shock capturing techniques is comparatively smaller than that dedicated to linear stabilization. As shown in the present paper, linear stabilization techniques can promote the Gibbs phenomenon and can transform a converging method into a non-converging one. We think that nonlinear stabilization methods (shock capturing, nonlinear viscosities, etc.) deserve more attention. These methods should be considered the workhorses that kill the Gibbs phenomenon and ensure convergence to the physically-correct weak solution; linear stabilization techniques should be seen as auxiliary tools whose job is to improve the convergence of an already converging nonlinear stabilization method. We refer to [START_REF] Burman | On nonlinear artificial viscosity, discrete maximum principle and hyperbolic conservation laws[END_REF] for a recent step in this direction.

Appendix A. Entropy viscosity. The purpose of this appendix is to briefly recall the principles of the so-called entropy viscosity method. The specificities of the method are unrelated to the content of this paper, but we nevertheless describe how it is implemented for the sake of completeness. Full descriptions of the technique can be found in [START_REF] Guermond | Entropy-based nonlinear viscosity for Fourier approximations of conservation laws[END_REF][START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF]. The time approximation of (2.6) is done using an explicit Runge-Kutta method defined by a Butcher tableau The algorithm is initialized by setting ν 0 h = 0 and u 0 h = u 0,h . The viscosity field is piecewise constant over the mesh K h . Assume that the current time is t n and let u n h be the approximation of u h (t n ). The next time increment is defined by (A.2) ∆t n = CFL min

K∈K h h K |f ′ (u n h | K )|
.

The approximation u n+1 h and the viscosity field ν n+1 h are evaluated at time t n+1 = t n + ∆t n as

(A.3) u n+1 h = u n h -∆t n (b 1 k 1 + b 2 k 2 + • • • + b r k r ), ν n+1 h = ν r ,
where the functions {k l } 1≤l≤r and the viscosity fields {ν l } 1≤l≤r are recursively computed as follows:

(A.4) The functions H 2 , . . . , H r are piecewise constant over the mesh K h and are defined as follows. First, we define the so-called entropy functional E ∈ C 1 (R; R). Second, we evaluate the entropy residual associated with two given states c h , d h ∈ X h as follows. Let δ be the time increment separating the states c h and d h . The entropy residual in the cell K at x ∈ K is defined as We now define the entropy residual associated with the mesh interfaces. Let F ∈ F i h . Then, for all x ∈ F , we set (A.6)

                                       w 1 = u n h , ν 1 = ν n h , Ω k 1 v dΩ =
J F (d h )(x) = (f ′ (d h (x))•n(x))[[∂ n E(d h (x))]].
The total entropy residual in the cell K is then defined to be (A.7) D K (c h , d h , δ) = max( R K (c h , d h , δ) L ∞ (K) , max

F ∈∂K∩F i h J F (d h ) L ∞ (F ) ).
The maximum local wave speed in the cell K is estimated by

(A.8) β K (d h ) = | max x∈K f ′ (d h (x))|.
Then, the value over K of the viscosity function H l , 2 ≤ l ≤ r, is defined as follows:

(A.9) H l (w 1 , . . . , w l )| K = min(c max β K (w l )h K , c ev h 2 K D K (w 1 , w l , c l ∆t n )) Note that each step of (A.4) requires solving a mass matrix linear system. We stress again that we do not recommend to use the lumped mass matrix instead of the consistent mass matrix since mass lumping induces large dispersion errors.

Figure 2 . 1 .

 21 Figure 2.1. One-dimensional linear transport. Results at t F = 1 without and with edge stabilization on five uniform meshes composed of 100, 2×100, . . . , 2 4 ×100 cells.

  Edge stabilized sol.

Figure 2 . 2 .

 22 Figure 2.2. Non-convex flux. Results at t F = 1 on a uniform mesh composed of 1000 cells and CFL = 0.01. Galerkin (left) and edge stabilized (right) solutions.

Figure 2 . 3 .

 23 Figure 2.3. Non-convex flux. Results at t F = 1 with entropy viscosity without edge stabilization (left) and with edge stabilization on a uniform mesh composed of 1000 cells (center and right).

Figure 2 . 4 .

 24 Figure 2.4. Left: One-dimensional Burgers equation, t F = 0.25, CFL = 0.025, uniform mesh composed of 200 cells; first-order viscous solution without (in green) and with edge stabilization (in red). Center and right: Nonconvex flux, t F = 1, CFL = 0.025, uniform meshes composed of 4, 000 (in red) and 10, 000 cells (in green); first-order viscous solution with edge stabilization (center: c ed = 1 in (2.9); right: c ed = 0.05 in (2.16)), the entropic solution is shown in dashed blue.

Figure 3 . 1 .KFigure 3 . 2 .

 3132 Figure 3.1. Definition of ∆ F (grey triangles)

  .

  3.

Figure 3 . 3 .

 33 Figure 3.3. The faces in thick red lines are in the set F ♯ h ; the triangles filled in dark grey are in Ω ∆♯ ; the triangles filled in dark or light grey are in Ω ♯ (and thus belong to K ♯ h ); unfilled triangles belong to K ♭ h .

Figure 4 . 1 .

 41 Figure 4.1. Nonconvex flux. Results at t F = 1 with entropy viscosity with weighted edge stabilization on five uniform meshes composed of 100, 200, 400, 800, and 1600 cells.

  Entropy visc., cmax = 1 2 , cev = 1.

  Ent. visc. + edge stab., cmax = 1 2 , cev = 1.

  Ent. visc. + weighted edge stab., cmax = 1 2 , cev = 1.

Figure 4 . 2 .

 42 Figure 4.2. Two-dimensional conservation equation with non-convex flux, t F = 1. Entropy viscosity solution without edge stabilization (left), with unweighted edge stabilization (center), and with weighted edge stabilization (right) on a uniform mesh composed of 118850 P 1 nodes, CFL = 0.025.

c

  r a r1 a r2 • • • a r,r-1 b 1 b 2 • • • b r-1 b r

Ω

  v∇•f (w 1 ) + b ev (ν 1 , w 1 , v), ∀v ∈ X h w 2 = u n h -∆t n a 21 k 1 , ν 2 = H 2 (w 1 , w 2 ), Ω k 2 v dΩ = Ω v∇•f (w 2 ) + b ev (ν 2 , w 2 , v), ∀v ∈ X h . . . w r = u n h -∆t n (a r1 k 1 + a r2 k 2 + • • • + a r,r-1 k r-1 ) ν r = H r (w 1 , w 2 , . . . , w r ), Ω k r v dΩ = Ω v∇•f (w r ) + b ev (ν r , w r , v), ∀v ∈ X h ,where b ev (ν, w, v) := K∈K h ν| K K ∇w•∇v dK.

(A. 5 )

 5 R K (c h , d h , δ)(x) = E(d h (x)) -E(c h (x)) δ + f ′ (d h (x))•∇(E(d h (x))).

Table 2 . 1

 21 One-dimensional linear transport. Maximum and minimum of entropy viscosity solution (left) and entropy viscosity + edge stabilized solution (right);P 1 finite elements, t F = 1.elements on various uniform meshes, and we evaluate the indicators e Max and e Min defined by (2.12).

		(a) Entropy viscosity		(b) Entropy viscosity + edge stabilization
	h	eMin	rate	eMax	rate	eMin	rate	eMax	rate
	2.500E-03	6.725E-03	-	6.715E-03	-	1.597E-02	-	1.597E-02	-
	1.250E-03	5.441E-03 0.306 5.434E-03 0.305	1.600E-02	-0.003 1.600E-02 -0.003
	6.250E-04	2.855E-03 0.930 2.854E-03 0.929	1.633E-02	-0.030 1.633E-02 -0.030
	3.125E-04	2.235E-03 0.353 2.235E-03 0.353	1.626E-02	0.006	1.626E-02	0.006
	1.563E-04	1.785E-03 0.324 1.785E-03 0.324	1.646E-02	-0.017 1.646E-02 -0.017

* This material is based upon work supported by the National Science Foundation grant DMS-0510650 and by the GNR MOMAS (CNRS/PACEN, ANDRA, BRGM, CEA, EDF, IRSN). Draft version, February 27, 2012

Remark 3.1. Theorem 3.1 holds in any space dimension under the assumption of locally quasiuniform meshes provided the parameter δ = sup h>0 max K∈K h max K ′ ∩K =∅ |1 -(h K ′ /h K ) 2 | is small enough, since the estimates η L ∞ (Ω) h k+1 and h -1 2 η L 2 (Ω) h k+ 1 2 hold, see [START_REF] Boman | Estimates for the L 2 -projection onto continuous finite element spaces in a weighted Lp-norm[END_REF]. We now prove a convergence result in three dimensions (with a slightly less stringent regularity assumption on the exact solution). The main difference with the result of Theorem 3.1 is that, at least for the lowest-order polynomials, the convergence result now requires that the parameter λ controlling the decrease of the weighting function α at infinity be not too large.

Theorem 3.2. Let u and u h be the solutions to (3.1) and to (3.2), respectively. Assume that u ∈ C 1 (0, t F ; H k+1 (Ω)) ∩ C 0 (0, t F ; W 1,∞ (Ω)). Assume d = 3 and quasi-uniform meshes. Then, the following holds for all t ∈ [0, t F ]: if k ≥ 3 and any λ,

and the right-hand side becomes h 2+ 2-λ 4+λ if k = 2 and λ < 2, and h 1+ 2-3λ 4+λ if k = 1 and λ < 2 3 . Proof. We proceed as in the proof of Theorem 3.1, up to the bound on T 2,2 (η, e) in step 3 of the previous proof. Here, we no longer use approximation properties in L ∞ (Ω) of the L 2 -orthogonal projection. Instead, owing to the Sobolev embedding, we observe that there holds η L p (Ω) h k with p = 6. Then,

and bounding meas(Ω ♯ ) as before yields

Collecting the above estimates and dropping higher-order terms, we arrive at

For k ≥ 4, we can take ǫ = 0 so that, with p = 6, 2k -1 + p-2 p (k -1 + ǫ) = 2k -1 + 2 3 (k -1) ≥ 2k + 1, yielding the desired estimate. In the case k ≤ 3, we use a bootstrap argument. Suppose that the error estimate (3.8) holds with bound h k+ρn . Then, the new bound for T 2,2 (η, e) is h 2k-1+ρn+ p-2 p (k-1+ǫn+ρn) , the bound on T 2,1 (η, e) remaining unchanged. The two terms are balanced by the choice

and this leads to the new error estimate with bound h k+ρn+1 where ρ n+1 = 1 2 (1λǫ n ). Hence, we obtain a recursive relation of the form ρ n+1 = a + bρ n . We want the fixed-point iteration to deliver a positive value. This, in turn, requires a > 0. Moreover, the desired

, which is larger than 1 2 only for k = 3. This completes the proof. Remark 3.2. Sharper results can be derived by using the advective derivative in the weighting function. Following [START_REF] Guermond | Subgrid stabilization of Galerkin approximations of linear contraction semi-groups of class C 0 in Hilbert spaces[END_REF], the proof of the error estimate then hinges on discrete inf-sup stability, assuming additional regularity in time of the exact solution. As a result, the desired h k+ 1 2 error estimate is achieved for any polynomial degree, any value for λ, and shape-regular meshes.

4. Numerical experiments. The objective of this section is to illustrate the efficiency of the weighted edge stabilization introduced and analyzed in §3.

4.1. One-dimensional tests. We report in this section one-dimensional tests over the periodic domain Ω := (0, 1) using P 1 finite elements. The time stepping is done with the SSP RK3 scheme with CFL = 0.2. The parameters for the first-order viscosity, entropy viscosity, and edge stabilization are c max = 1 2 , c ev = 1 2 , and c ed = 0.05, respectively. 4.1.1. Smooth transport. We consider the one-dimensional transport equation ∂ t u + ∂ x u = 0 with initial data u 0 (x) = sin(2πx) to compare the convergence properties of the linear edge stabilization with those of the weighted edge stabilization. The computations are performed up to t F = 1 on various meshes. The meshes are non-uniform to avoid super-convergence effects. The size of each cell is chosen randomly with deformation factor equal to 3, that is, the size ratio between two neighboring cells is at most 3.

The L 1 -and L 2 -norms of the error at t F = 1 are reported in Table 4.1 for seven meshes. This test shows that the convergence properties of the linear edge stabilization and those of the weighted edge stabilization are identical for smooth solutions. 4.1.2. Non-smooth transport. We now test the convergence properties of the entropy viscosity method with the weighted edge stabilization. We consider the one-dimensional transport equation with non-smooth initial data (2.11) and periodic boundary conditions. The computations are performed with P 1 finite elements up to t F = 1 on various meshes. The L 1 -and L 2 -norms of the error are reported in Table 4.2 for seven uniform (left table) and seven non-uniform meshes (right table). The convergence orders in the L 1 -and L 2 -norm seem to be 3 4 and 3 8 , respectively. Considering that the initial data is in 2 ) for smooth solutions. We now verify that the weighted edge stabilization does preserve the weakened maximum principle (see (

Two-dimensional tests.

We illustrate the method on two-dimensional problems.

Smooth transport.

Let us consider the domain Ω = {x ∈ R 2 , x < 1}, the rotating velocity field β(x) = 2π(-y, x), where x = (x, y), and the two-dimensional transport problem

with a = 0.3 and r 0 = (0.4, 0). Note that ∇•β = 0 and β•n| ∂Ω = 0.

We perform tests on triangular isoparametric finite elements. The meshes are quasi-uniform and based on a Delaunay triangulation technique [START_REF] Rebay | Efficient unstructured mesh generation by means of Delaunay triangulation and Bowyer-Watson algorithm[END_REF]. The time stepping uses the RK4 scheme and the CFL used in all the tests is CFL=0.25. The parameters for the first-order viscosity, entropy viscosity, and edge stabilization are c max = 1 4k , c ev = 0.1, and c ed = 0.025, respectively. Convergence tests on the above smooth problem, not reported here, have shown that, when used alone, the weighted edge stabilization performs as well as the original edge stabilization technique both for P 1 and P 2 finite elements. When combined with the weighted edge stabilization, the converge rate of the entropy viscosity method with P 1 finite elements, which is already second-order in the L 2 -norm, is not increased. The story is a little different for higher-order finite elements. The entropy viscosity alone has been observed to be second-order for P 1 finite elements and to be of order (k + ǫ(k)) for higher-order polynomial degrees; see [START_REF] Guermond | Entropy viscosity method for nonlinear conservation laws[END_REF]. Augmenting the entropy viscosity method with weighted edge stabilization improves the convergence order. We show in Table 4.4 convergence tests on the above smooth transport problem (4.1) with P 2 finite elements. We compare the convergence rate of the entropy viscosity method and the entropy viscosity method with weighted edge stabilization. Adding the weighted edge stabilization clearly improves the convergence rate of the entropy viscosity method. Some super-convergence effect is observed (the errors are estimated by using a Gaussian quadrature exact on P 5 polynomials). 4.4 Linear transport with smooth solution (4.1). L 1 -and L 2 -norms of error with the entropy viscosity method (left) and the entropy viscosity with weighted edge stabilization (right), P 2 finite elements, t F = 1.

4.2.2.

Non-smooth transport. We solve again the linear transport problem (4.1), but this time we consider the following non-smooth initial data u(x, 0) = 1 if xr 0 < a and u(x, 0) = 0 otherwise. The meshes and the time stepping are the same as before.

We report in Table 4.5 convergence tests in the L 1 -and L 2 -norms for various meshes using the entropy viscosity method with weighted edge stabilization. The results for P 1 elements are shown in the left table and those for P 2 elements are shown in the right table. The convergence orders in the L 1 -and L 2 -norms seem to be 3 4 and 3 8 for P 1 elements and 4 5 and 2 5 for P 2 elements, respectively. As above, these rates are compatible with the estimates 2 ) for smooth solutions. We now verify the weakened maximum principle by computing the indicators e Max and e Min defined by (2.12) for each mesh. The results are reported in Table 4 ∂ t u + ∇•f (u) = 0, f (u) = (sin u, cos u), u(x, y, 0) = 3.5π, x 2 + y 2 < 1; 0.25π, otherwise.

The local velocity is f ′ (u) = (cos u,sin u). This is a Cauchy problem in R 2 . This test was proposed in [START_REF] Kurganov | Adaptive semidiscrete central-upwind schemes for nonconvex hyperbolic conservation laws[END_REF] and is challenging to many high-order numerical schemes because the solution has a two-dimensional composite wave structure; see [START_REF] Kurganov | Adaptive semidiscrete central-upwind schemes for nonconvex hyperbolic conservation laws[END_REF] for details.

We perform computations on a mesh composed of triangular finite elements. The mesh family is quasi-uniform. The solution is computed at t F = 1. The time stepping is done with the standard RK4 scheme with CFL = 0.025. Three computations are done with the entropy viscosity method (with coefficients c max = 1 2 , c ev = 1): one without edge stabilization, one with unweighted edge stabilization (with c ed = 1), and one with weighted edge stabilization (with c ed = 1 and λ = 2). The graph of the P 1 solution without edge stabilization is shown in Figure 4.2(a). The solution shows the expected rotating composite wave structure composed of a shock and an expansion. This solution is an accurate approximation of the entropic solution. The graph of the solution with unweighted edge stabilization is shown in Figure 4.2(b). We observe that this solution exhibits a non-physical


Charles Bourassin-Bouchet

Sébastien De Rossi
email: sebastien.derossi@institutoptique.fr

Jihuan Wang

Evgueni Meltchakov

Angelo Giglia

Nicola Mahne

Stefano Nannarone

Franck Delmotte

Shaping of single-cycle sub-50-attosecond pulses with multilayer mirrors

Keywords: numbers: 42, 65, Re 42, 65, Ky 07, 85, Fv

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

Introduction

The perfect synchronization of a broadband emission in the extreme ultraviolet (XUV) range is the key to generate ultrashort pulses at the attosecond time scale. Such a radiation can nowadays be produced thanks to a process named high order harmonic generation (HHG) [1,2], resulting from the nonlinear interaction of an infrared (IR) femtosecond field with a gaseous [3,4] or solid target [5]. However, the spectral components of such attosecond pulses were proved to be linearly desynchronized [6], due to a positive group delay dispersion (GDD). This phenomenon commonly called the atto-chirp is intrinsic to HHG, and inevitably broadens pulses. To compensate for this temporal chirp, it has been proposed either to optimize the HHG conditions [7], or to add on the way an optical component inspired by what already existed for visible or IR femtosecond pulses. Such a post-compression can be perform thanks to chirp multilayer mirrors [8,9,10,11], grating compressors [12], or the natural dispersion of metallic filters [13] or gaseous media [14]. Some of these solutions have been successfully implemented experimentally and allowed one to get isolated 80 as pulses [15] and trains of 63 as pulses [14]. The phenomenon preventing from obtaining even shorter pulses appears to be the width of the spectral range on which the chirp compensation is performed.

However, in the case of visible femtosecond pulses, a solution called the timedomain optimization, has been developed to control the synchronization of a broadband radiation using chirped multilayer mirrors. Indeed, instead of targeting a wanted spectral behaviour, i. e. a wanted GDD, during the design of the multilayer structure, one directly optimizes the temporal profile of the reflected pulse. This improves the compression efficiency of the mirrors, which allowed one to get single-cycle sub-5 f s pulses [16,17]. It was proposed to transpose this approach into the XUV domain [18] to get 100 as pulses.

In this article, we report on the theoretical and experimental study of time domain optimized multilayer mirrors to obtain sub-50 as pulses. To reach such a short duration, the dispersion has to be controlled over a spectral range larger than 100 eV . This bandwidth is isolated from the rest of the spectrum by an aluminum-beryllium filter. Then the resulting pulse is shaped by the optimized multilayer mirror. We confirmed the robustness of the pulse profiles for reasonable experimental conditions. We then fabricated and characterized the mirrors in amplitude and in phase on synchrotron radiation. The spectral phases of the mirrors were reconstructed over a 120 eV bandwidth using the total electron yield standing wave method [19]. Finally, we deduced the attosecond reflected electric field from these measurements.

Conditions for the optimization

The attosecond source

In order to design a multilayer mirror suitable for realistic conditions, we chose typical generation parameters for the simulations of the spectrum S source (ω) and the spectral phase ϕ source (ω) of the attosecond source. To be more specific, the wavelength of the driving femtosecond laser pulse equals 800 nm, and its intensity I IR is chosen to be 5.8 × 10 14 W.cm -2 , which is a typical value [6,14,15]. The XUV radiation is generated in neon. Moreover, the well-known cut-off law [20] predicts that the cut-off energy of an HHG spectrum equals I p +3.17 U p , with I p the ionization potential of the gas and U p the ponderomotive energy of the IR pulse. In these conditions, the cut-off energy is equal to 132 eV , the GDD of ϕ source in the plateau region equals 3600 as 2 .rad -1 , and becomes zero in the cut-off region [21]. Moreover, when plotted in logarithm scale, experimental spectra often appear to decrease linearly in the plateau region [6,10], and to rapidly fall down in the cut-off region. Thus, we used this simplified shape for S source (ω) (inset of figure 1 (b) as a basis for the optimization. Using equation (1), we deduced that the attosecond pulses I source (t) obtained in these conditions have a duration equal to 132 as and a typical profile [6] composed of a main peak followed by rapidly decreasing [13], a 300 nm beryllium filter [14], a 300 nm zirconium filter [15], and a 150 nm palladium filter [11]. Transmissions are calculated using the CXRO data base above 30 eV . Below, the Palick indices are used. This explains the imperfect connection of the transmission of the Al filter at 30 eV . Inset: The transmission of the Al-Be filter and the normalized XUV spectrum (red curve) are plotted in logarithm scale. The change of slope on the spectrum stands for the cut-off region.

bounces, see figure 2 (a).

I source (t) = +∞ -∞ S source (ω) 1/2 exp (iϕ source (ω)) • exp (iωt) dω 2 (1) = F T S 1/2 source • exp (iϕ source ) 2 2.

The filter

The filter is a key component since it isolates the wanted spectral range from the rest of the spectrum. Especially, by absorbing lower energies, it allows one to suppress the IR beam used for the generation of the XUV radiation. Furthermore, the characterization of an attosecond pulse is usually based on the measurement of photoelectron spectra coming from the ionization of a rare gas by the XUV pulse [3,4]. As a consequence, only photons carrying an energy greater than the ionization potential of the gas, say 15.76 eV in argon, will be represented in the electron spectra. Thus the use of a filter allows one to select a spectral range detectable in typical attosecond experiments.

Based on these considerations, we first considered typical filters, as depicted in figure 1 (b). As is well-known, the larger the spectrum, the shorter the Fourier-limited pulse. But it rapidly appeared that none of these filters had a bandwidth large enough to attain sub-50 as durations with the chosen spectrum. Nevertheless, if combining together the 17 -73 eV bandwidth of an aluminum filter with the 40 -112 eV bandpass of a beryllium filter, one can select a very large range. However, adding the two independent filters in the beam will transmit the radiation only in the range where the two bandwidths overlap. But if their thicknesses are weak enough, one can get a significant transmission on the complete bandwidth. Consequently, we combined in a single filter a 70 nm thick Be filter with a 40 nm thick Al filter. Since the obtained hybrid filter is only made of a few layers and is used in transmission, the interferential effects in the stack are negligible. So the complex spectral response of the filter t f ilter is equivalent to the product of the Al and Be filters taken independently. In order to ensure that the spectral response will be the same in the two possible ways of the filter, we made it symmetrical by encapsulating the 70 nm thick Be layer between two 20 nm Al layers. We also added two 5 nm aluminum oxides Al 2 O 3 layers on each side of the filter as in a realistic case, see the inset of figure 1 (a). Thus, the chosen filter has a total thickness greater than 100 nm and is easily manufacturable [22]. Figure 1 (b) shows that its transmission simulated using the CXRO and Palick indices is comparable to other typical filters over the 17 -112 eV range. Moreover, its transmission falls down below 17 eV , as depicted in the inset of figure 1 (b), and equals 1.9 × 10 -5 and 2.5 × 10 -7 at 800 nm and 11 eV respectively. This ensures a good rejection of the IR beam and of the lowest energies. Finally, figure 2 (b) shows that the pulse I in (t) (see equation ( 2)) transmitted by the filter is 103 as long, whereas the duration of the source pulse was equal to 132 as. This compression comes likely from the natural dispersion of aluminum and berrylium, which slightly compensates for the atto-chirp. Moreover, the duration of the Fourier limited pulse obtainable after the filter equals 47 as, which means that the bandpass of the filter is large enough to get sub-50 as pulses. However, these durations could be attained only if the phase dispersion is compensated by a properly designed multilayer mirror.

I in (t) = F T S 1/2 source • exp (iϕ source ) • t f ilter 2 = F T S 1/2 in • exp (iϕ in ) 2 (2)

The multilayer mirrors

Multilayer mirrors for the XUV range are typically based on a periodic stack of molybdenum (Mo) and silicon (Si) layers. But being given that the radiation incoming on the mirror has complex spectral properties, the mirror spectral response optimizing the pulse shape and duration is not intuitive. Consequently, the best multilayer stack could be arbitrarily complex, that is to say not necessarily periodic. Thus we decided to optimize totally aperiodic structures. Moreover, in order to increase the reflectivity at low energy [23] and to minimize the interdiffusion of layers, we added a boron carbide (B 4 C) layer between each Si and Mo layers. So the structure is composed of Mo-B 4 C-Si-B 4 C elementary stacks repeated several times. Due to the significant absorption of these materials in the selected spectral range, we expected that the total number of Transform limited pulse (dashed line). For a), the spectrum below the 17 eV has been put to zero to compare the pulses before and after the filter on the same spectral range. c) Multilayer structure of the mirror ML 1 . (d-f) Final pulse profile after the reflection off ML 1 , ML 2 and ML 3 , respectively. For each case, the FWHM duration, the relative fluence η, and the normalized intensity of the secondary peak are reported.

layers should not exceed a few tens of layers. Since the materials and the order of the layers were already chosen, the optimization process consisted only in finding their thicknesses. However, in order to obtain realistic multilayer stacks, we imposed the minimal thickness of a layer to be 0.8 nm. Finally, we added on top of the stack a 2 nm thick Si layer, to allow the phase measurement (see section 4.), and a 1 nm thick SiO 2 layer to take into account the unavoidable oxidation process. We chose an interfacial roughness equal to 0.3 nm in the entire stack, which is usually obtained in our deposition machine. The last parameter to consider was the angle of incidence of the XUV beam on the structure. Due to the total reflection phenomenon occurring for grazing angles, the larger the incidence angle, the greater the reflectivity, but the less comprehensive the control of the phase of the mirror. Thus, we made a trade-off and chose to optimize the multilayer structures for a 22.5 • incidence angle.

The temporal optimization algorithm

The optimization procedure of the thicknesses was based on a usual simulated annealing algorithm [24]. During the numerical process, the algorithm calculated the complex reflectivity of the multilayer r mir (ω) as a function of the photon energy hω. Then, knowing r mir and the incoming spectrum S in and spectral phase ϕ in , the algorithm calculated the electric field of the reflected pulse I out (t) using equation ( 3):

I out (t) = F T S 1/2 in • exp (iϕ in ) • r mir 2 = F T S 1/2 out • exp (iϕ out ) 2 (3)
As it is the case for every optimization algorithm, the simulated annealing finds the solution by minimizing an objective function. For the time-domain optimization, it had to take into account the relevant parameters of the wanted pulse. We chose to optimize the full-width at half maximum (FWHM) duration τ , which is the duration systematically given when measuring an isolated attosecond pulse [4,15,25] or an attosecond pulse train [3,14]. However, the FWHM only describes the width of the main peak in the pulse profile, and ignores the rest of the pulse. Thus, we minimized the relative highness I 2nd of the second highest peak in the pulse profile in order to only get a main burst in the envelope. Finally, to make sure that the optimized pulses would carry a reasonably high energy, we chose the normalized fluence η (see equation ( 4)) of the pulse as the third parameter of the merit function, as shown in equation ( 5):

η = S out (ω) dω/ S in (ω) dω (4)
OF = w τ • (τ -τ target ) 2 + w I 2nd • (I 2nd -I target 2nd ) 2 + w η • (η -η target ) 2 (5)
Where the X target indicates the value targeted for the parameter X during the optimization, and w X is the corresponding weight of the parameter X.

Optimization results and robustness

Results

We optimized three multilayer structures referred to as ML 1,2,3 . The stacks are made of about twenty layers, the structure of the mirror ML 1 being reported on figure 2 (c) as a typical case. The pulse profiles obtained after reflection off each mirror are depicted in figure 2 (d-f). It appears that they are all composed of a main peak carrying the major part of the energy in the middle of a weakly intense 1 f s temporal pedestal. This typical pulse shape is a direct consequence of the set of parameters chosen for the objective function, since neither τ , nor η, nor I 2nd is able to limit the temporal extension of this pedestal. Moreover, the relative fluence of the optimized pulses is greater than 1% in the three cases, which remains high enough to perform a pulse characterization [15].

Regarding the pulse duration and the second peak intensity, two approaches were chosen. In the case of ML 1 , we imposed that I 2nd had to be as low as possible, and andc) the intensity of the IR source generating the attosecond pulses were chosen as variable parameters. In each case, the continuous pulse evolution with respect to the chosen parameter is depicted (central panel), the pulse profile being detailed for three values of the parameter (upper panels). The evolution of the FWHM and the intensity of the secondary peak I 2nd (lower panels) is also depicted.

relaxed the condition on the FWHM which had to be lower than 60 as. Concerning ML 2 and ML 3 , the trade-off was made on the second peak intensity which had to be lower than 20%. The targeted FWHM was the duration of the Fourier limited pulse after the filter (see the dashed line on figure 2 (b) ), namely 47 as. We see on figure 2 (e-f) that the latter duration was attained for both mirrors and that I 2nd is equal to 16% and 10% for ML 2 and ML 3 respectively. Finally, if comparing the pulse generated by the source and the three pulse profiles obtained after the mirrors, we see that the FWHM has been divided by almost three.

Robustness to environmental conditions

However, these theoretical pulse profiles were optimized assuming given initial conditions, such as the shape of the initial XUV spectrum, or the angle of incidence of the pulse on the mirror. To make sure that the obtained mirrors will perform a pulse compression even though the environmental conditions slightly change, we simulated the evolution of the pulse profiles with respect to three parameters. The robustness of the three mirrors evolving the same way, we represented only ML 3 in order not to overload figure 3. First of all, in order to take into account the impact of a realistic fabrication process of the mirrors, we considered a deviation of the thicknesses of the deposited layers from the expected ones. To be more specific, we randomly changed the thicknesses of each layer of the targeted multilayer structure of ML 3 , the amplitude of the perturbation being given by the peak-to-valley (PV) thickness variation . In order to get a reliable tendency, we randomly drew 100 sets of thicknesses for each value of . Among the 100 obtained pulses, we plotted the one with the worst τ and I 2nd values in the upper and central parts of figure 3 (a), and we calculated the mean of τ and I 2nd , (see the lower panel of figure 3 (a) ). The results highlight that should not exceed a few angströms in order to get the wanted pulse profile. Hopefully, the typical precision of our deposition machine is close to 0.3 nm PV on the thickness of each individual layer, which lies in this tolerance range.

As a second parameter, we considered the incidence angle of the XUV beam on the mirror. The mirrors being optimized for an angle of 22.5 • , we plot in figure 3 (b) the evolution of the pulse profile on the [0 • ; 45 • ] range. The simulations show that the tolerance range of the mirror is remarkably large. Indeed, a compression of the pulse below 55 as can be achieved in normal incidence as well as for a 30 • incidence angle while keeping I 2nd below 20%. Moreover, the compression becomes optimal between 17.5 • and 25 • , the duration being below 50 as and I 2nd being lower than 15%. Above 30 • , the total reflection phenomenon becomes problematic. Indeed, the grazer the beam, the less it penetrates in the multilayer structure, the less effective the shaping of the pulse. Such a range of possible angles of incidence makes a pulse compression possible in various experimental setups.

The last variable parameter was chosen to be the intensity of the IR beam generating the XUV pulses. This quantity determines the cut-off energy of the HHG process, i. e. the extension of the XUV spectrum, and the value of the GDD responsible for the attochirp [6]. To explore a wide range around the initial intensity equal to 5.8×10 14 W.cm -2 , we changed I IR from 1.5 × 10 14 W.cm -2 to 15 × 10 14 W.cm -2 . Consequently, the cutoff energy shifted from 50 eV to 306 eV , and the GDD decreased from 14000 as 2 .rad -1 down to 1400 as 2 .rad -1 . According to the simulations, the optimal range goes from 5 × 10 14 W.cm -2 to 8.5 × 10 14 W.cm -2 . In these conditions, the FWHM remains below 50 as and I 2nd does not exceed 20%. In an experimental set up, the typical precision when setting the IR intensity is much better than 1 × 10 14 W.cm -2 [6], so it is reasonable to think that one can easily work in the optimal range of compression.

These simulations showed that the pulse compression was possible in realistic experimental conditions. Consequently, the optimized multilayer structures were relevant and could be manufactured.

Fabrication and characterization of the mirrors

Fabrication

The multilayer stacks were deposited using a magnetron sputtering machine equipped with four cathodes [23]. Successive layers were deposited by scanning the substrate above the target at a distance of 10 cm. During the process, we used a 0.27 P a argon pressure in the deposition chamber. The plasma discharges was established with a RF power of 150 W for B 4 C and Si targets and a DC current of 0.06 A for the Mo target. The multilayer structure was controlled by means of a grazing-incidence Cu-Kα reflectometer. Finally, the experimental performances of the deposited mirrors have been measured at the wavelengths of interest on a synchrotron beamline.

Principle of the characterization

Being given that the incident spectrum is very broad, the mirrors spectral response had to be characterized over a very large spectral range, namely from the cut-on of the Al-Be filter at 17 eV up to the cut-off energy of the initial spectrum at 132 eV . In this range, the complete characterization of the mirrors implied to have an access both to the squared modulus R mir and the phase ϕ mir of the complex reflectivity r mir . If the measurement of R mir is now common, the determination of the phase is not straightforward. However, a possible method is based on the measurement of the Total Electron Yield, i. e. the current coming from the sample irradiated with synchrotron monochromatic light. In first approximation, this current is proportional to the standing wave coming from the interference between the amplitudes of the incident beam A in (ω) and of the reflected beam A out (ω) on top of the stack [19,26]. In S polarization, the measured photocurrent i mir (ω) is given by equation ( 6):

i mir (ω) = C (ω) • |A in (ω) + A out (ω)| 2 = C (ω) • |A in | 2 • 1 + R mir + 2R 1/2 mir • cos (ϕ mir ) (6)
In this equation, C (ω) is a parameter depending on the material from which the electrons are emitted, i. e. the material of the surface layer. To overcome this unknown, one needs to measure a reference photocurrent i ref (ω). The latter has to come from a reference sample identical to the top layer of the multilayer stack, and for which C (ω) is supposed to be the same. Since the surface layer of the ML mirrors was chosen to be a silicon layer, the reference sample had to be composed of a simple silicon layer. Thus, if evaluating J (ω) corresponding to the ratio i mir (ω) i ref (ω) , one can remove the influence of C (ω) as shown by equation (7).

J (ω) = i mir (ω) i ref (ω) = 1 + R mir + 2R 1/2 mir • cos (ϕ mir ) 1 + R ref + 2R 1/2 ref • cos (ϕ ref ) (7)
But, if the reflectivities of the mirror R mir (ω) and of the reference sample R ref (ω) are easily measurable, the reference phase ϕ ref (ω) remains unknown. However, if rearranging equation (7) into equation (8), one can notice that the lower the ratio R ref R mir , hereafter referred to as ρ, the less significant the role of ϕ ref in the reconstructed phase ϕ mir .

cos (ϕ mir (ω)) = J • 1 + R ref 2R 1/2 mir + ρ • cos (ϕ ref ) - 1 + R mir 2R 1/2 mir (8)
Moreover, cos (ϕ ref ) being always included in the [-1; 1] interval, one can limit cos (ϕ mir ) by two functions L ± using inequality (9). Consequently, measuring J (ω), R mir (ω) and R ref (ω) allows one to reconstruct the phase of the multilayer mirror with an accuracy depending on the quantity ρ.

L -(ω) ≤ cos (ϕ mir (ω)) ≤ L + (ω) (9)
Where

L ± equals J 1 + R ref ± 2R 1/2 ref -1 + R mir • 2R 1/2 mir -1
.

Results of the characterization

According to the previous analysis, the photoelectrons are supposed to come only from the surface of the measured material. Consequently, i ref (ω) should not depend on the thickness of the silicon reference sample. To verify this hypothesis, we manufactured a 10 nm thick and a 20 nm thick silicon reference sample, and measured their reflectivities and photocurrents on the BEAR beamline on the Elettra synchrotron, see figure 4 (a) and (b). As expected, the photocurrents of the two reference samples are extremely similar. As for their reflectivities, they are very low in both cases, and are modulated with an oscillation, the period of which is determined by the thickness of the silicon layer. After the fine characterization of the reference samples, it was possible to use them to extract the phase of the ML mirrors. Consequently, we measured the reflectivity and the photocurrent of the three mirrors, and used the phase extraction process previously described. The reference phase ϕ ref being unknown, we plotted in the case of ML 1 the range of the possible cos (ϕ mir ) using equation ( 9), as depicted in figure 4 (c). When ρ is low, roughly on the 40 -100 eV spectral range, the shaded area delimited by L -and L + strongly narrows, allowing one to estimate the phase of the mirror with a good precision. Moreover, due to the modulations in the reflectivities of the reference samples, R ref and thus ρ become periodically very low. On the opposite, for the lowest and highest energies, R ref becomes of the same order of magnitude than R mir , see figure 4 (a). Consequently, ρ becomes significant leading to a bad accuracy of the phase measurement. Applying the arccos function to the data in figure 4 (c) allowed one to get the phase ϕ mir of the mirror. The reconstructed phase had then to be unwrapped by changing its sign and by adding a π offset every time the phase reaches 0 or π. Finally, the two experimental unwrapped phases reconstructed from the two reference samples were averaged, the final phase being depicted on figure 5 (d-f) for the three mirrors along with their reflectivities in figure 5 (a-c).

Concerning ML 1 and ML 2 , the agreement between the experimental phase and the simulated one is very good, the latter being almost always included in the shaded area. Moreover, the global tendency of the phase as well as its most rapid variations are well reproduced in the 40 -100 eV range. The same agreement can be observed for their reflectivity too, except around 40 eV , as shown in figure 5 (a) and (b). This is most likely due to a deviation of the actual indices of refraction in this range from the tabulated ones used for the optimization. Regarding ML 3 , the agreement appears to be less good.

If the global shape of the reflectivity is more or less respected, the two main spikes around 80 eV are clearly shifted of a few electron-volts with respect to the simulated reflectivity. Moreover, the same shift affects the cosine of the phase extracted from the photocurrents measurements, see the inset in figure 5 (f), and the unwrapped phase.

After studying the multilayer mirror by the use of a Cu-Kα reflectometer, it appears that this shift comes from the formation of a new chemical compound at each interface between a B4C layer and a Mo layer, when these two layers are too thin. Finally, it should be noticed that the simulated and measured spectral responses of these mirrors are not intuitive. Until now, when designing mirrors for the compression of attosecond pulses, most of the targeted phase behaviors were a constant and negative GDD, that is to say a convex parabolic phase profile [9,10,11,19,26]. With the temporal optimization, the obtained phase profiles are not parabolic anymore. More precisely, the phases of the three mirrors clearly tend to be convex, but they are modulated by fluctuations around this global curvature. However, the latter do not prevent from obtaining an effective compression, the duration of the final pulse being guaranteed by the time domain optimization. This observation has already been made for chirped dielectric mirrors for the compression of IR femtosecond pulses [16], and appears to be relevant for the compression of attosecond pulses too.

In order to completely characterise the mirrors, the ultimate step is to reconstruct the reflected pulses from the experimental responses of the mirrors. To do so, one needs to consider a given phase evolution in the range of the possible phases, i. e. to choose a curve in the shaded area in figure 5 (d-f). Thus, we chose the case where cos (ϕ ref ) equals zero as an average case, see the inset of figure 5 (d). The resulting pulse profiles are depicted in figure 6. In the three cases, the reconstructed pulses are very close to the expected ones. Especially, the durations obtained after ML 2 and ML 3 are below 50 as while keeping the second peak below 20 %.

After reflection off one of the mirrors, the attosecond pulses become single-cycle, as shown on figure 7. In an experiment involving such pulses, the response of the medium interacting with the pulse will not be averaged over many optical cycles. The medium will rather interact directly with the electric field of the pulse. And this interaction will strongly depend on the carrier-envelope phase (CEP) of the pulses, that is the phase-shift between the centre of the pulse and the maximum of the electric field. Moreover, a theoretical study [25] highlighted that the CEP of attosecond pulses was robust to reasonable variations of the generation conditions. Consequently, attosecond pulses could be generated without significant shot to shot variations of their CEP state. However, if a post-compression of these attosecond pulses is performed with multilayer mirrors, a modification of their CEP can be induced during the reflection off the mirror. Thus, it becomes important to estimate experimentally the impact on the attosecond CEP. Furthermore, the carrier-envelope phase cannot be measured by current attosecond pulse characterization techniques, such as the RABITT technique [3,6,10] or the CRAB technique [11,27]. Indeed, these methods reconstruct the pulse spectral phase up to a constant, the information about the CEP being encoded in this very constant. However, the modification of the CEP induced by the mirror is accessible through the photocurrents technique, since the latter allows one to reconstruct the absolute spectral phase of the mirror. Consequently, we made a study based on the experimental results to estimate the variation of the CEP after the reflection off the mirrors.

Control of the Carrier-Envelope Phase of attosecond pulses

Rigorously, one can define the CEP ϕ CE of a pulse as the opposite of the temporal phase φ (t) taken at the arrival time t 0 of the pulse, that is when the maximum intensity is reached :

ϕ CE = -φ (t 0 ) (10)
Thus the first step when determining the CEP of a pulse is to identify its arrival time t 0 .

Time spent by the pulse in the multilayer stack

Consider a fictitious detector placed at a given position along the pulse propagation axis. One can define the arrival time t 0 of the pulse as the moment when the pulse reaches its maximum intensity on the detector. If putting a multilayer mirror on the way, the variation of the arrival time ∆t 0 will correspond to the time that the pulse spent in the multilayer structure, see figure 7. As shown in figure 7 (a), the duration of the reflection measured in the case of ML 1 equals 245 ± 3 as, which is very close to its theoretical value of 243 as. Regarding the other mirrors, the pulse should spend more time in the multilayer stack since ∆t 0 was estimated to be equal to 327 ± 2 as and to 317 ± 4 as for ML 2 and ML 3 respectively.

Carrier-Envelope Phase offset induced by the mirror

But the reflection of the attosecond pulse on a mirror can also induce a change of its CEP. This phase shift ∆ϕ CE can be evaluated on figure 7. Indeed, in this figure, the CEP of the incident pulse was chosen to be zero, so that estimating the CEP of the reflected pulse is equivalent to measuring ∆ϕ CE . Especially, in the case of ML 1 , the measured CEP shift is equal to 2.2 ± 0.3 rad, the theoretical value being equal to 2.0 rad. As for ML 2 and ML 3 , we estimated ∆ϕ CE to be equal to -2.2 ± 0.1 rad and to -1.9 ± 0.4 rad respectively. This highlights some similarities in the behaviour of ML 2 and ML 3 . Indeed, for these two mirrors, the duration of the reflection process is slightly longer than in the case of ML 1 , and the CEP shift that they induce appears to be negative, as opposed to a positive shift for ML 1 . This measurement of the offset of the CEP due to the reflection on a mirror has to be seen as a proof of concept which will have to be improved. Moreover, even though the influence of the mirror on the CEP can be measured, it has not been taken into account in the optimization process. Thus we can imagine modifying the objective function to optimize a mirror compressing attosecond pulses and giving them a wanted CEP state.

Conclusion

By using the time domain approach, we optimized a broadband filter and three multilayer mirrors to compress attosecond pulses down to theoretical durations between 59 as and 47 as FWHM. Moreover, a theoretical study highlighted that the tolerance of the mirrors to the variations of the experimental conditions were good. In particular, the mirrors can ensure an efficient pulse compression in a very large range of angles of incidence. We then manufactured the mirrors and characterized them on a synchrotron beamline both in reflectivity and phase. Especially, the measurement of the latter was performed thanks to the photocurrent technique over a 120 eV large spectral range. The very good overall agreement with simulated phases allowed us to estimate experimentally the time spent by an attosecond pulse in the multilayer mirrors during the reflection process. We also measured the impact of the mirrors on the attosecond carrier-enveloppe phase. This opens the route toward attosecond electric field shapers based on multilayer mirrors.

Figure 1 .

 1 Photon energy (eV)

Figure 2 .

 2 Figure 2. Theoretical pulse profiles before (a) and after (b) the filter, and after (d-f ) reflection off one of the three multilayer mirrors optimized in the time-domain. a) (resp. b) Pulse profiles (continuous line) generated by the attosecond source (resp. transmitted by the filter) and the associated FourierTransform limited pulse (dashed line). For a), the spectrum below the 17 eV has been put to zero to compare the pulses before and after the filter on the same spectral range. c) Multilayer structure of the mirror ML 1 . (d-f) Final pulse profile after the reflection off ML 1 , ML 2 and ML 3 , respectively. For each case, the FWHM duration, the relative fluence η, and the normalized intensity of the secondary peak are reported.

Figure 3 .

 3 Figure 3. Evolution of the pulse profile after reflection off ML 3 with respect to environmental parameters. a) The PV thickness error of the deposited layers, b) the angle of incidence of the XUV beam on the mirror,and c) the intensity of the IR source generating the attosecond pulses were chosen as variable parameters. In each case, the continuous pulse evolution with respect to the chosen parameter is depicted (central panel), the pulse profile being detailed for three values of the parameter (upper panels). The evolution of the FWHM and the intensity of the secondary peak I 2nd (lower panels) is also depicted.

Figure 4 .

 4 Photon energy (eV) Photocurrent (arb. u.)

Figure 5 .

 5 Photon energy (eV) 60 100 120

Figure 6 .

 6 Figure 6. Reconstructed pulse profiles based on the measured mirrors properties. (a-c) The shaded lines correspond to the pulse profiles deduced from the measured reflectivities and from the phases obtained if assuming cos (ϕ ref ) equal to 0. The theoretical profiles (dashed line) are also represented. In each case, the duration FWHM, the second peak intensity and the normalized fluence η are depicted along with the theoretical values.

Figure 7 .

 7 Figure 7. Experimental and theoretical influence of the mirrors on the reflected attosecond electric field. For ML 1 (a), ML 2 (b) and ML 3 (c), the incident pulse (purple shaded line) and the reconstructed reflected pulses (red, yellow and blue shaded lines) are represented. The continuous oscillating lines stand for the experimentally reconstructed electric field after reflection off the mirrors, the dashed lines corresponding to the simulated electric field. As shown in the inset in a), the CEP of the incident pulse was assumed to be 0. In each case, the experimental and theoretical values of the duration of the reflection of the pulse ∆t 0 , and of the CEP shift ∆ϕ CE are depicted.