Olivier Lafitte

Olof Runborg

ERROR ESTIMATES FOR GAUSSIAN BEAMS AT A FOLD CAUSTIC

In this work we show an error estimate for a first order Gaussian beam at a fold caustic, approximating time-harmonic waves governed by the Helmholtz equation. For the caustic that we study the exact solution can be constructed using Airy functions and there are explicit formulae for the Gaussian beam parameters. Via precise comparisons we show that the pointwise error on the caustic is of the order O(k -5/6 ) where k is the wave number in Helmholtz.

Introduction

Gaussian beam superpositions is a high frequency asymptotic approximation for solutions of wave equations [START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF]. It is used in numerical methods to simulate waves in the high frequency regime. Unlike standard geometrical optics, the Gaussian beam approximation does not break down at caustics, which is one of its main advantages.

In this paper we consider error estimates for the approximation in terms of the wave number k > 0. Error estimates for Gaussian beams are known in a number of settings. See for instance [START_REF] Liu | Error estimates for Gaussian beam superpositions[END_REF][START_REF] Liu | Sobolev and max norm error estimates for Gaussian beam superpositions[END_REF] and the references therein. The main result is that, in L 2 and Sobolev norms, the relative error of first order beams decays as O(k -1/2 ), independently of dimension and regardless of the presence of caustics. This has been shown for general strictly hyperbolic partial differential equations and the Schrödinger equation [START_REF] Liu | Error estimates for Gaussian beam superpositions[END_REF][START_REF] Liu | Sobolev and max norm error estimates for Gaussian beam superpositions[END_REF] as well as the Helmholtz equation [START_REF] Liu | Gaussian beam method for the Helmholtz equation[END_REF]. The better rate O(k -1 ) is typically observed in numerical computations and has been shown in L 2 for the Schrödinger equation [START_REF] Zheng | Optimal error estimates for first-order Gaussian beam approximations to the Schrödinger equation[END_REF], and also in L ∞ for the Schrödinger and the acoustic wave equation on sets strictly away from caustics [START_REF] Liu | Sobolev and max norm error estimates for Gaussian beam superpositions[END_REF]. Similar estimates have been also been derived for higher order beams. For p-th order beams the rates are O(k -p/2 ) and O(k -p/2 ) respectively. There are, however, no precise, pointwise, error estimates for the solution at a caustic. In particular, for first order beams it has not been shown that this error vanishes as k → ∞, although there is ample numerical evidence to this effect.

The purpose of this paper is to show such an error estimate for a typical fold caustic in two dimensions. More precisely, we consider the Helmholtz equation [START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF] ∆u + k 2 (1 -x)u = 0.

We assume there is an incident wave u inc from x = -∞ making an angle Θ ∈ (0, π/2) with the x-axis. Moreover, at x = 0 it has the amplitude envelope A(y), so that u inc (0, y) = A(y)e ik sin Θ .

This wave will generate a fold caustic at the line x = x c where

x c = cos 2 (Θ).
Figure 1 shows a representative solution. In Section 3 we make this physical situation precise and derive an exact solution using Airy functions on R 2 . We subsequently study the solution in the region 0 ≤ x ≤ x c and compare it at x = x c to an approximation using Gaussian beams, denoted by u GB . (Note that 0 < 1 -x ≤ 1 in this region; we do not make comparisons elsewhere, as the equation then no longer models the physical situation.) The main result is the following theorem. Theorem 1.1. Suppose A is a Schwartz class function, A ∈ S(R), and that 0 < Θ0 ≤ Θ ≤ Θ1 < π/2. Let u be the exact solution of (1) as defined in Section 3 and u GB a first order Gaussian beam approximation detailed in Section 4. Then there is a constant C independent of k such that

||u GB (x c , •) -u(x c , •)|| L ∞ ≤ Ck -5/6 .
This result hence confirms that first order Gaussian beams do converge pointwise at the caustic. Moreover, since the solution itself grows as O(k 1/6 ) at this caustic, the relative error is O(k -1 ), the same as away from the caustic. We conjecture that this will be the case also for more general caustics.

The paper is organized as follows. In Section 2 notations are established and some preliminary results are discussed. In Section 3 the exact solution is defined and a formula for it is derived. In Section 4, the corresponding Gaussian beam approximation is introduced. Sections 5, 6 and 7 contain estimates of the Gaussian beam parameters, the phase, various oscillatory integrals, as well as the exact solution and the Gaussian beam approximation. In Section 8 the proof of the main result in Theorem 1.1 is carried out. Finally, in Section 9 some properties of Airy functions are presented.

Preliminaries

In the analysis we use a k-scaled Fourier transform and indicate it with a hat mark on the function, [START_REF] Liu | Error estimates for Gaussian beam superpositions[END_REF] f (η) := F k (f )(η) := √ kF(f )(kη) = k 2π f (x)e -ikxη dx.

The corresponding scaled inversion formula reads

f (x) = k 2π f (η)e ikxη dη =: F -1 k ( f )(x) = √ kF -1 ( f )(kx).
We also have

f x (η) = ikη f (η), ||f || L ∞ ≤ k 2π || f || L 1 , F k (f * g)(η) = 2π k f (η)ĝ(η)
We will frequently make use of a smooth, even, cut-off function which we denote ψ ∈ C ∞ c (R). It is defined as

ψ(x) =      1, |x| ≤ 1, 0, |x| ≥ 2, ∈ (0, 1), 1 < |x| < 2, ψ(-x) = ψ(x).
This is used to divide integrals into subdomains and to regularize the Fourier transform of functions in S , the space of tempered distributions. For example, if f is in L ∞ (R), but not in L 1 (R), the definition (2) must be interpreted in distributional sense. We then let ψ t := ψ(x/t) and consider instead the Fourier transform of the compactly supported function f ψ t , which is well-defined by (2) for all t > 0. The following Lemma shows that the limit as t → ∞ gives us the Fourier transform in S . Lemma 2.1. Let f ∈ L ∞ and set ψ t = ψ((b 0 x + b 1 )/t) for any fixed real numbers b 0 = 0 and b 1 . Then, with F k as defined above,

lim t→∞ F k (f ψ t ) = F k (f ),
in S . Moreover, if g ∈ S then, again in S ,

lim t→∞ F k ((f * g)ψ t ) = lim t→∞ 2π k F k (f ψ t )F k (g) = 2π k F k (f )F k (g).
The short proof is found in the Appendix. In particular, if f is defined pointwise, the Lemma shows that [START_REF] Liu | Sobolev and max norm error estimates for Gaussian beam superpositions[END_REF] f (η) = lim t→∞ k 2π ψ(x/t)f (x)e -ikxη dx.

We also introduce some notation that will prove useful later on in the paper. We let (4) ξ 0 = cos Θ, η 0 = sin Θ, so that (ξ 0 , η 0 ) T is the unit vector pointing in the propagation direction of the incident wave. Following Theorem 1.1 we will assume, throughout the paper, that 0 < Θ0 ≤ Θ ≤ Θ1 < π/2. This translates to bounds on ξ 0 and η 0 of the form

(5) 0 < ξ0 ≤ ξ 0 ≤ ξ1 < 1, 0 < η0 ≤ η 0 ≤ η1 < 1.
for some ξj and ηj . Moreover, we let

(6) δ = x c -x = ξ 2 0 -
x, be the distance to the caustic. Finally, we introduce the polynomial q, which is related to the geometrical spreading of the rays, [START_REF] Ludwig | Uniform asymptotic expansions at a caustic[END_REF] q(s) = 1 + 2is -s 2 β, β = 1 + 2iξ 0 .

It will be used frequently in the analysis.

Expression of the exact Helmholtz equation solution

In this section we define an exact solution to the Helmholtz equation for the physical setup described in the introduction. Using a property of the Airy function we deduce a decomposition of the solution into forward and backward going waves.

We consider a solution u to (1) which is a tempered distribution on R 2 , i.e. u ∈ S (R 2 ). The solution then has a k-scaled Fourier transform in y which we denote û(x, η). Upon Fourier transforming also [START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF] in y, we obtain an ODE for û(x, η),

(8) ûxx + k 2 (1 -x -η 2 )û = 0.
The only tempered distribution solution to this ODE is given by

û(x, η) = a(η)Ai(k 2/3 (x + η 2 -1)),
where Ai is the Airy function of the first kind, and a(η) is a function to be determined. This solution is thus a C ∞ bounded solution. The Airy function in this expression contains waves going both forward and backward. In the sequel, we will choose the function a(η) as k 1/6 P (k, η) defined in [START_REF] Lafitte | The kernel of the Neumann operator for a strictly diffractive analytic problem[END_REF]. To arrive at this choice, we first note that when α = exp(iπ/3), it holds for all z that [12, Eq. 9.2.14], [START_REF] Lebeau | Regularite gevrey 3 pour la diffraction[END_REF] Ai(z) = αAi (-αz) + ᾱAi (-ᾱz) .

This follows since Ai(-αz) and Ai(-ᾱz) solve the same ODE as Ai(z) given that α 3 = -1. We then introduce the scaled variables

ζ(x) = k 2 3 (x + η 2 -1), ζ + (x) = -αζ, ζ -(x) = -ᾱζ, such that û(x, η) = a(η)Ai(ζ(x)) = a(η)αAi(ζ + (x)) + a(η)ᾱAi(ζ -(x)).
To further understand this decomposition, we note that the asymptotics of the Airy function (in the angular sector |arg(z

)| < 2π/3) is Ai(z) 1 2 √ π z -1 4 exp - 2 3 z 3 2
.

Therefore, upon defining the phase φ(x)

= 2 3 (1 -η 2 -x) 2/3 we have, when x < 1 -η 2 , Ai(ζ + (x)) Ck -1/6 e -ikφ(x) (1 -η 2 -x) 1 4 , Ai(ζ -(x)) Ck -1/6 e +ikφ(x) (1 -η 2 -x) 1 4
.

Thus, the phases of the two expressions, and their gradients, have opposite signs, meaning that the two terms in (9) represent a decomposition into forward and backward going waves in the region x < 1 -η 2 . More precisely, the solution is fully known when a(η) is given, and one assumes that this solution of ∆u + k 2 (1 -x)u = 0 for all x (including x < 0 where the velocity grows) is the sum of an incoming and an outgoing wave of the form a(η)αAi(ζ + (x)) and a(η)ᾱAi (ζ -(x)), respectively. This leads us to define

T 0 (x, η) = Ai (ζ(x)) Ai (ζ(0)) , T + (x, η) = Ai (ζ + (x)) Ai (ζ + (0)) , T -(x, η) = Ai (ζ -(x)) Ai (ζ -(0)) .
which are three particular solutions of ( 8), since each term in (9) solve the ODE. These solutions are normalized such that they equal one for x = 0. The solutions T + and T -represent forward and backward going waves. Among the three solutions, only T 0 is bounded, since Ai(z) is bounded, while T ± are not even in S , since Ai(αz) includes also the unbounded second kind Airy function Bi(z), cf. (72). We are looking for the solution û(x, η) = T 0 (x, η)û(0, η).

We do not know û(0, η), just that the incoming part of û(0, η) represents the incident plane wave. We therefore write it as a sum of an incoming and an outgoing part û(0, η) = û+ (0, η) + û-(0, η),

where we define û+ (0, η) as the k-scaled Fourier transform in y of the incoming wave with amplitude A and direction Θ (recall η 0 = sin Θ),

u + (0, y) = u inc (0, y) = A(y)e ikη0y .
We then want to find û-(0, η) such that

T 0 (x, η)û(0, η) = T + (x, η)û + (0, η) + T -(x, η)û -(0, η).
To achieve this, it is necessary and sufficient that the values of the functions and their derivatives agree at x = 0, since both sides satisfy the same second order ODE. This gives us the linear relations

û(0, η) = û+ (0, η) + û-(0, η), T 0 (0, η)û(0, η) = T + (0, η)û + (0, η) + T -(0, η)û -(0, η), from which we can deduce û-(0, η) = T + (0, η) -T 0 (0, η) T 0 (0, η) -T -(0, η) û+ (0, η) =: T (η)û + (0, η). It follows that û(x, η) = T 0 (x, η)û(0, η) = T 0 (x, η)(û + (0, η) + û-(0, η)) = T 0 (x, η)(1 + T (η))û + (0, η).
The decomposition is not valid at the roots of Ai(ζ(0)) = Ai(k

2 3 (η 2 -1)
). Another form is available, which is valid at all points. It is given in the following Lemma. Lemma 3.1. One can express T (η) as follows

T (η) = -α Ai(ζ -(0)) Ai(ζ + (0)) , 1 + T (η) = ᾱ Ai(ζ(0)) Ai(ζ + (0)) , Proof.
To simplify notation we write Ai +,-,0 (x) := Ai(ζ +,-,0 (x)). Then, using (9),

T ± -T 0 = Ai ± (x) Ai ± (0) - Ai 0 (x) Ai 0 (0) = Ai ± (x)Ai 0 (0) -Ai ± (0)Ai 0 (x) Ai ± (0)Ai 0 (x) = Ai ± (x)(αAi + (0) + ᾱAi -(0)) -Ai ± (0)(αAi + (x) + ᾱAi -(x)) Ai ± (0)Ai 0 (x) = ±α ∓1 Ai -(0)Ai + (x) -Ai + (0)Ai -(x) Ai ± (0)Ai 0 (x) .
Hence,

T + -T 0 = (T --T 0 ) Ai -(0) Ai + (0) (-ᾱ 2 ) = (T --T 0 ) Ai -(0) Ai + (0) α.
The first identity follows upon differentating this expression with respect to x. The second identity is then given by another application of [START_REF] Lebeau | Regularite gevrey 3 pour la diffraction[END_REF].

It follows now that

T 0 (x, η)(1 + T (η)) = ᾱ Ai(ζ(x)) Ai(ζ(0) Ai(ζ(0)) Ai(ζ + (0)) = ᾱ Ai(ζ(x)) Ai(ζ + (0)) .
Since we know û+ (0, η) we can thus express the full solution as

û(x, η 0 + η) = ᾱ Ai(k 2 3 (x -X)) Ai(αk 2 3 X) û+ (0, η 0 + η),
where we defined

X(η) = 1 -(η 0 + η) 2 .
In this expression, one notices that the denominator never vanish because all the roots of the Airy function are on the negative real axis.

Finally, since u + (0, y) = A(y)e ikη0y ⇒ û+ (0, η) = Â(η -η 0 ), we get û(x,

η 0 + η) = ᾱ Ai(k 2 3 (x -X)) Ai(αk 2 3 X) Â(η).
We write this as û(x, η 0 + η) = v(η, x, k) Â(η), where [START_REF] Lafitte | The kernel of the Neumann operator for a strictly diffractive analytic problem[END_REF] v(η, x, k) = k 1/6 P (k, η)Ai(k 2 3 (x -X)), P (k, η) = ᾱk -1/6 Ai(αk 2/3 X) .

Construnction of the Gaussian beam approximation

In this section we derive expressions for a first order Gaussian beam approximation to the solution of [START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF]. A Gaussian beam is a high frequency asymptotic solution to the Helmholtz equation. To model a general solution of (1), superpositions of Gaussian beams are used. We give the general form of a Gaussian beam and their superposition in R 2 below. The derivations of the expressions can be found in [START_REF] Liu | Gaussian beam method for the Helmholtz equation[END_REF].

The Helmholtz equation with a general index of refraction n(x) reads [START_REF] Wasow | Linear turning point theory[END_REF] ∆u + k 2 n(x) 2 u = 0.

When n(x) is real, the equation models wave propagation, but it has a well-defined solution also when n(x) is imaginary. However, Gaussian beams can only be defined for real n(x).

A first order Gaussian beam for [START_REF] Wasow | Linear turning point theory[END_REF] has the form

(12) v b (x) = a(s)e ik(S(s)+(x-γ(s))•p(s)+ 1 2 (x-γ(s)) T M (s)(x-γ(s)) , s = s * (x),
where a(s) ∈ C is the amplitude, S(s) ∈ R the reference phase, p(s) ∈ R 2 the phase gradient and M (s) ∈ C 2×2 the phase Hessian. Moreover, γ(s) ∈ R 2 is the central ray, which agrees with the rays of geometrical optics. An example of a Gaussian beam is shown in Figure 2. In (12) the parameter s depends on the point of evaluation x. Normally one takes the s-value for the point on the central ray that is closest to x, as indicated in Figure 2. However, in the analysis below we make a simpler choice. By a result in [START_REF] Ralston | Gaussian beams and the propagation of singularities[END_REF] the Hessian M will always have a positive definite imaginary part. The solution v b will therefore be a "fattened" version of the central ray, with a Gaussian profile normal to the ray with a width determined by M . The s-dependent parameters in the Gaussian beam are all given by ODEs [START_REF] Liu | Gaussian beam method for the Helmholtz equation[END_REF], as follows

dγ ds = 2p, dp ds = ∇n 2 (γ), dM ds = D 2 (∇n 2 )(γ) -2M 2 , ( 13
)
dS ds = 2n 2 (γ), da ds = -tr(M )a.
The initial data for γ and p is given by the starting point (x 0 , y 0 ) and direction (ξ 0 , η 0 ) of the beam,

γ(0) = x 0 y 0 , p(0) = ξ 0 η 0 .
In order to form an admissible Gaussian beam, M (0) must always satisfy

(14) M (0) T = M (0), M (0)γ (0) = p (0), a T (ImM (0))a > 0, when a ⊥ γ (0).
The choice of M (0) and the precise form of the incoming wave finally determine the initial data for S and a. We come back to this issue below.

To build more general solutions we use superpositions of Gaussian beams. We assume that the incoming wave is known along a curve Γ in R 2 , which we can parameterize with the parameter z, so that Γ(z) = (x 0 (z), y 0 (z)). For each point on Γ we launch one Gaussian beam in the direction θ 0 (z) of the wave at that point. The parameters of the beams then also depend on z and we write γ = γ(s; z), p = p(s; z), etc. This gives the beams v b = v beam (x; z), from which we finally construct the Gaussian beam superposition

(15) u GB (x) = k 2π v beam (x; z)dz.
See [START_REF] Liu | Gaussian beam method for the Helmholtz equation[END_REF] for more details. In a numerical scheme the z-variable is discretized and for each discrete value, the ODEs (13) are solved with a numerical ODE method. The superposition (15) is subsequently computed using numerical quadrature.

4.1. Expressions for Gaussian beam parameters. In (1) we have the index of refraction n(x) = √ 1 -x. The ODEs (13) can be solved explicitly and we get analytic formulae for all parameters in the Gaussian beam. To show this we let p = (ξ, η) and γ = (x, y) and also recall that p(0) = (ξ 0 , η 0 ), where ξ 2 0 + η 2 0 = 1. Since we let the beam start at x 0 = 0, we set γ(0) = (0, y 0 ). With the particular choice of n it follows from (13) that ξ (s) = -1 and η (s) = 0. Hence,

ξ(s) = ξ 0 -s, η(s) = η 0 .
For the positions, we get x (s) = 2ξ(s) = 2ξ 0 -2s and y (s) = 2η(s) = 2η 0 . By also using the initial data we obtain x(s) = 2sξ 0 -s 2 , y(s) = y 0 + 2sη 0 .

The caustic x = x c is located at the point where the ray turns back, i.e. where x (s) = 0, which gives s = ξ 0 and x c = x(ξ 0 ) = ξ 2 0 . Note that all the rays are confined to the region x ≤ x c < 1, where the index of refraction is real-valued. The fact that n(x) is complex-valued for x > 1 therefore does not affect the Gaussian beams.

We also need to compute the coefficients corresponding to the phase S, the second derivative of the phase M and the amplitude a. We have

dS ds = 2n 2 (x(s)) = 2(1 -x(s)) = 2[1 -2sξ 0 + s 2 ], S(0) = S 0 ,
so the phase is a third order polynomial,

S(s) = S 0 + 2s -2s 2 ξ 0 + 2 3 s 3 .
For M we have the Riccati equation

dM ds = D 2 (n 2 )(x(s)) -2M (s) 2 = -2M (s) 2 , M (0) = M 0 ,
with the solution M (s) = (I + 2sM 0 ) -1 M 0 .

The matrix M 0 must satisfy the conditions in (14). We pick

M 0 = Q + iP, P = η 2 0 -η 0 ξ 0 -η 0 ξ 0 ξ 2 0 , Q = 1 2 -ξ 0 -η 0 -η 0 ξ 0 .
Note that P, Q are symmetric, P is the orthogonal projection on p ⊥ 0 , and 2Qp 0 = -e 1 = p (0). Moreover,

(p ⊥ 0 ) T Qp ⊥ 0 = ξ 0 > 0. Next, one checks that M (s) = (I + 2sM 0 ) -1 M 0 = 1 q(s) (I + 2s(iI -M 0 ))M 0 , where det(I + 2sM 0 ) = 1 + 2sTrM 0 + (2s) 2 detM 0 = 1 + 2is -s 2 β = q(s), β = 1 + 2iξ 0 .
We note that q is related to the geometrical spreading of the beams. Further manipulations, using the facts that P 2 = P , 4Q 2 = I and 2(P Q + QP -Q) = ξ 0 I reveals that M (s) can be written simply as We let m ij be the elements of M and deduce that

M (s) = 1 q(s) M 0 - 1 2 sβI .
m 11 (s) = - 1 2q(s) (ξ 0 -2iη 2 0 + sβ) = - 1 2q(s) (-2i + (ξ 0 + s)β), ( 16
)
m 22 (s) = 1 2q(s) (ξ 0 + 2iξ 2 0 -sβ) = 1 2q(s) (ξ 0 -s)β. (17) Finally, for a, da ds = -tr(M (s))a, a(0) = a 0 .
We note that if λ 0 and λ 1 are eigenvalues of M 0 , then q(s)

= det(I + 2sM 0 ) = (1 + 2sλ 0 )(1 + 2sλ 1 ) and tr(M (s)) = λ 0 1 + 2sλ 0 + λ 1 1 + 2sλ 1 = 1 2 d ds log q(s).
It follows that (18) a(s) = a 0 q(s) .

The last thing needed to make the expression (12) for the Gaussian beam well defined, is to decide which s-value to use for a given x, i.e. the function s * (x). As mentioned above, this is normally taken to be the s-value for the point on the central ray that is closest to x. Here, however, to simplify, we just take the s-value for the point of the curve that has the same y-coordinate; see Figure 3. With x = (x, y) this leads to

s * (x) = s * (y) = y -y 0 2η 0 . Then x -γ(s * (x)) = (x -x(s * ), 0) T and (12) becomes v b (x, y) = a (s * ) exp ik S(s * ) + (x -x(s * ))ξ(s * ) + 1 2 m 11 (s * )(x -x(s * )) 2 ,
with x(s), ξ(s), S(s), m 11 (s), a(s) and s * (y) given above.

Remark 4.1. The value of q(s) = 1 + 2si -s 2 β crosses the negative real axis when s = 1/ξ 0 . To find a better branch cut for the square root in the expression (18) for A(s) we note first that the equation

q(s) β = 0 ⇔ ( β)( q(s)) = ( β)( q(s)) ⇔ 2ξ 0 (1 -s 2 ) = 2s(1 -ξ 0 s)
has the unique solution s = ξ 0 . Therefore, the equation q(s) = tβ with t ∈ R only has a solution if t = q(ξ 0 )/β = (1 -ξ 2 0 ) > 0. Hence, q(s) never crosses the line {-tβ : t ≥ 0}, which we therefore use as branch cut. This guarantees a smooth dependence of the Gaussian beam on s for all s ≥ 0. It can be written as

√ zβ * / √ β * if √
• is the usual square root with branch cut along the negative real axis.

4.2.

Gaussian beam superposition. We will now prepare the superposition. The initial curve Γ is simply the y-axis so that Γ(z) = γ(0; z) = (x 0 (z), y 0 (z)) = (0, z). By assumption the incoming wave propagates in the same direction (ξ 0 , η 0 ) for all z. Moreover, the same initial data for M is used for all z. This means that x, ξ, η, M and the geometrical spreading parameter q are all independent of z. Only y, S and a depend on z. We obtain from (12),

v beam (x, y; z) = a(s; z)e ik(S(s;z)+(x-γ(s;z))•p(s)+ 1 2 (x-γ(s;z)) T M (s)(x-γ(s;z)) , s = s * (x).
To derive the initial data a 0 (z), S 0 (z) for a(s; z) and S(s; z) we consider the trace of v beam and u GB on x = 0. To get explicit formuale we also let s * (x) be the s-value for the point on the curve that has the same x-value, i.e. s * (0, y) = 0. That gives

v beam (0, y; z) = a 0 (z)e ik[S0(z)+((0,y)-γ(0;z))•p 0 + 1 2 ((0,y)-γ(0;z)) T M (0)((0,y)-γ(0;z))] = a 0 (z)e ik[S0(z)+(y-z)η0+ 1 2 (y-z) 2 m22(0)] .
and for u GB ,

u GB (0, y) = e ikyη0 k 2π a 0 (z)e ik[S0(z)-zη0]+ 1 2 ik(y-z) 2 m22(0) dz = a 0 (y)e ikyη0 k 2π e ik[S0(z+y)-(z+y)η0]+ 1 2 ikz 2 m22(0) dz + O(k -1 ).
To match this with the incoming wave on x = 0, i.e. u inc (0, y) = A(y) exp(ikη 0 y), we take

S 0 (z) = η 0 z, and
a 0 (y) = A(y) k 2π e 1 2 ikz 2 m22(0) dz -1 = A(y) -im 22 (0).
Thus the expressions for the Gaussian beam coefficients are

x(s) = 2sξ 0 -s 2 , (19a) y(s; z) = z + 2sη 0 , (19b) ξ(s) = ξ 0 -s, (19c) η(s) = η 0 , (19d) S(s; z) = η 0 z + 2s -2s 2 ξ 0 + 2 3 s 3 , (19e) m 11 (s) = 2i -(ξ 0 + s)β 2q(s) , (19f) a(s; z) = A(z) q(s) -im 22 (0). (19g)
This gives us the simplified expression for v beam , (20)

v beam (x, y; z) = a (s * ; z) e ik(S(s * ;z)+(x-x(s * ))ξ(s * )+ 1 2 m11(s * )(x-x(s * )) 2 ) , s * (y; z) = y -z 2η 0 , which, together with (19), define u GB via (21) u GB (x, y) = k 2π v beam (x, y; z)dz.
We will continue now to simplify (20) and compute the k-scaled Fourier transform of u GB in y. Since x -x(s * ) = x -2s * ξ 0 -s * 2 and δ = ξ 2 0 -x by ( 6) we have

v beam (x, y; z) = A(z)f (y -z)e ikη0y , f (2η 0 s * ) = -im 22 (0) e ikSg(s * ;δ) q(s * )
where

S g (s * ; δ) := S(s * ; z) + (x -x(s * ))ξ(s * ) + 1 2 m 11 (s * )(x -x(s * )) 2 -η 0 y = η 0 (z -y) + 2s * -2s * 2 ξ 0 + 2 3 s * 3 + (x -2s * ξ 0 + s * 2 )(ξ 0 -s * ) + 1 2 m 11 (s * )(x -2s * ξ 0 + s * 2 ) 2 = 2ξ 2 0 s * -2s * 2 ξ 0 + 2 3 s * 3 + ((s * -ξ 0 ) 2 -δ)(ξ 0 -s * ) + 1 2 m 11 (s * )((s * -ξ 0 ) 2 -δ) 2 = 2 3 ξ 3 0 + 2 3 (s * -ξ 0 ) 3 + ((s * -ξ 0 ) 2 -δ)(ξ 0 -s * ) + 1 2 m 11 (s * )((s * -ξ 0 ) 2 -δ) 2 .
Then we can write

u GB (x, y) = k 2π (A * f )(y)e ikη0y .
Since A ∈ S the beam v is always integrable in z, so that u GB in ( 21) is well-defined. However, there is no guarantee that u GB (x, •) is in L 1 ; in general it is not. When we compute its Fourier transform we therefore use Lemma 2.1 and (3). By Lemma 5.1 and 5.2 below, q is bounded away from zero and m 11 is strictly positive. Hence, f ∈ L ∞ . Therefore,

ûGB (x, η 0 + •) = lim t→∞ F k (ψ t u GB (x, •) exp(-ikη 0 •)) = lim t→∞ k 2π F k (ψ t (A * f )) = lim t→∞ F k (f ψ t )F k (A),
where we can choose ψ t (w

) := ψ((w/2η 0 -ξ 0 )/t), that is, b 0 = 1/2η 0 = 0 and b 1 = -ξ 0 in Lemma 2.1.
We then compute

F k (f ψ t )(δ, η) = k 2π ψ t (w)f (w)e -ikηw dw = k 2π ψ t (2η 0 s * )f (2η 0 s * )e -2ikηη0s * 2η 0 ds * = -im 22 (0) k 2π ψ t (2η 0 s * ) e ikSg(s * ;δ) q(s * ) e -2ikηη0s * 2η 0 ds * = -im 22 (0) k 2π ψ(θ/t) e ikSg(θ+ξ0;δ) q(θ + ξ 0 ) e -2ikηη0(θ+ξ0) 2η 0 dθ
Here we made the change of variables w = 2η 0 s * and θ = s * -ξ 0 . Moreover,

S g (θ + ξ 0 ; δ) = 2 3 ξ 3 0 + 2 3 θ 3 -(θ 2 -δ)θ + 1 2 m 11 (θ + ξ 0 )(θ 2 -δ) 2 = 2 3 ξ 3 0 - 1 3 θ 3 + δθ + 1 2 m 11 (θ + ξ 0 )(θ 2 -δ) 2 = 2 3 ξ 3 0 + φ g (δ, θ, η) + 2η 0 ηθ,
where we have introduced the Gaussian beam phase φ g as,

(22) φ g (δ, θ, η) = - 1 3 θ 3 + θ(δ -2η 0 η) + 1 2 m 11 (ξ 0 + θ)(θ 2 -δ) 2 .
Since m 22 (0) = ξ 0 β/2 = by (17) and q(ξ 0 ) = βη 2 0 , this finally gives

F k (f ψ t )(δ, η) = -im 22 (0) k 2π ψ(θ/t) e ik( 2 3 ξ 3 0 +φg(δ,θ,η)-2ηη0ξ0) q(θ + ξ 0 ) 2η 0 dθ = 2η 0 -im 22 (0) k 2π e ik( 2 3 ξ 3 0 -2ηη0ξ0) ψ(θ/t) e ikφg(δ,θ,η) q(θ + ξ 0 ) dθ = -2iξ 0 q(ξ 0 ) k 2π e ik( 2 3 ξ 3 0 -2ηη0ξ0) ψ(θ/t) e ikφg(δ,θ,η) q(θ + ξ 0 ) dθ =: k 1/6 P GB (k, η)I t (η, δ, k), where (23) P GB (k, η) = 2(πξ 0 ) 1/2 e -iπ/4 e ik( 2 3 ξ 3 0 -2ηη0ξ0) ,

and

(24)

I t (η, δ, k) = k 1/3 q(ξ 0 ) 2π ψ(θ/t) e ikφg(δ,θ,η) q(θ + ξ 0 ) dθ. Then (25) ûGB (x, η 0 + η) = vGB (η, x, k) Â(η), with vGB (η, x, k) = lim t→∞ k 1/6 P GB (k, η)I t (η, δ, k).

Properties of the amplitude and phases

In this section we collect a series of estimates that we will need for the final proof of the magnitude of the Gaussian beam error.

Geometrical spreading.

Here we show some properties of the q-polynomial in ( 7) that relates to the geometrical spreading, repeated here for convenience,

q(s) = 1 + 2is -βs 2 , β = 1 + 2iξ 0 .
We have Lemma 5.1. There are positive constants q 0 and q 1 , independent of ξ 0 and θ ∈ R, such that

0 < q 0 (1 + θ 2 ) ≤ |q(ξ 0 + θ)| ≤ q 1 (1 + θ 2 ).
Furthermore, there are constants b n independent of θ and ξ 0 such that

d n dθ n 1 q(ξ 0 + θ) ≤ b n 1 + |θ| n+1 .
Proof. We first show that q has no real root for the considered values of ξ 0 . Suppose therefore that q has a real root s = r. Then the real and imaginary parts of q(r) = 0 reads

1 -r 2 = 0, 2r -2ξ 0 r 2 = 0.
The only solutions to this system are r = ξ 0 = ±1, which are both ruled out by [START_REF] Zheng | Optimal error estimates for first-order Gaussian beam approximations to the Schrödinger equation[END_REF]. Let q(θ) := |q(ξ 0 + θ)|/(1 + θ 2 ), which is then continuous and non-zero for all θ. For large θ it is bounded from below and above since lim θ→±∞ q(θ) = |β|. In fact, there is a a constant θ 0 such that

|q(θ) -|β|| ≤ |β| 2 , |θ| > θ 0 ,
uniformly in ξ 0 , because of the bound (5). We can then take q 0 and q 1 as

q 0 = min inf |θ|≤θ0 q(θ), |β| 2 , q 1 = max sup |θ|≤θ0 q(θ), 3|β| 2 .
The stated bound then follows.

For the second statement, we observe that there exists a sequence of polynomials p n of degree n such that

d n dθ n 1 q(θ) = p n (θ) q(θ) n+1/2 ,
given by the recursion

p n+1 (θ) = p n (θ)q(θ) -(n + 1/2)p n (θ)q (θ), thanks to d dθ p n (θ) q(θ) n+1/2 = p n (θ)q(θ) -(n + 1/2)p n (θ)q (θ) q(θ) n+3/2 .
Then, by the lower bound on |q| and ( 5), there are constants C n and b n depending on n but independent of θ and ξ 0 such that

d n dθ n 1 q(ξ 0 + θ) ≤ |p n (ξ 0 + θ)| q n+1/2 0 (1 + θ 2 ) n+1/2 ≤ C n 1 + |ξ 0 + θ| n 1 + |θ| 2n+1 ≤ C n 1 + 2 n-1 (|ξ 0 | n + |θ| n ) 1 + |θ| 2n+1 ≤ b n 1 1 + |θ| n+1 .
This shows the lemma. 5.2. Phase. In this section we define the two phase functions φ a and φ g that turn up in our analysis and show a few properties of them. The first one, φ a , is defined as

(26) φ a (δ, θ, η) = - 1 3 θ 3 + θ(δ -2η 0 η).
This has a close connection to the Airy function and we call it the Airy phase. Indeed, exp(iθ 3 /3) is the Fourier transform of Ai(z) in S (R) and therefore, using the regularization of Lemma 2.1,

lim s→∞ k 1/3 2π ψ(θ/s)e -ikφa(δ,θ,η) dθ = Ai(k 2/3 (2η 0 η -δ)).
The second phase is the Gaussian beam phase (22) derived in the previous section. It can be written as a sum of φ a and a complex correction, by ( 16),

(27) φ g (δ, θ, η) = φ a (δ, θ, η) + 1 2 m 11 (ξ 0 + θ)(θ 2 -δ) 2 , m 11 (s) = 2i -(ξ 0 + s)β 2q(s) .
We start by looking at the m 11 part of the Gaussian beam phase.

Lemma 5.2. For m 11 it holds that

lim s→±∞ m 11 (s)s = lim s→±∞ -m 11 (s)s 2 = 1 2 , (28)
m 11 (s) = η 2 0 |q(s)| 2 , ( 29
)
d n m 11 (ξ 0 + θ) dθ n ≤ d n 1 + |θ| n+1 , ( 30
)
|m 11 (ξ 0 + θ)| ≥ D 0 1 + |θ| , ( 31
)
for some constants d n , D 0 independent of θ and ξ 0 .

Proof. From (27) we have

sm 11 (s) = s 2i -(ξ 0 + s)β 2q(s) = (2i -ξ 0 β)s -βs 2 2(1 + 2i -βs 2 ) = (-2i + ξ 0 β)s -1 + β -2(1 + 2i)s -1 + 2β → 1 2 ,
showing the first limit in (28). Moreover,

m 11 (s) = - 1 2 β(1 + 2is -s 2 β) + (2i -2sβ)(2i -(ξ 0 + s)β) q(s) 2 = - 1 2 β -4 -2iβξ 0 + 2sβ ξ 0 β -2i + s 2 β 2 q(s) 2
which similarly implies the second limit in (28). For (29) we have by ( 16),

m 11 (s) = (2i -(ξ 0 + s)β)(1 -2is -s 2 β * ) 2|q(s)| 2 = (ξ 0 + s)(2s -2s 2 ξ 0 ) + 2(η 2 0 -sξ 0 )(1 -s 2 ) 2|q(s)| 2 = η 2 0 + s (ξ 0 + s)(1 -sξ 0 ) -η 2 0 s -ξ 0 (1 -s 2 ) |q(s)| 2 = η 2 0 |q(s)| 2 > 0.
To show (31) with D 0 = min(1/2, 1 -ξ1 )/ √ 2q 1 we observe that,

|m 11 (ξ 0 + θ)| = |2i -(θ + 2ξ 0 )(1 + 2iξ 0 )| |q(θ + ξ 0 )| ≥ |θ + 2ξ 0 | + 2|1 -(θ + 2ξ 0 )ξ 0 | √ 2q 1 (1 + θ 2 ) ≥ |θ + 2ξ 0 | + 1 -|(θ + 2ξ 0 )ξ 0 | √ 2q 1 (1 + θ 2 ) = 1 + (1 -ξ 0 )|θ + 2ξ 0 | √ 2q 1 (1 + θ 2 ) ≥ 1 -2(1 -ξ 0 )ξ 0 + (1 -ξ1 )|θ| √ 2q 1 (1 + θ 2 ) ≥ D 0 1 + |θ| ,
where we used Lemma 5.1 as well as the facts that

√ 2|z| ≥ | z| + | z| and 1 -2x(1 -x) ≥ 1/2 for 0 ≤ x ≤ 1.
For (30) we note that if r k and p are any polynomials of degrees k and , then

d dθ r k p k = r k p -r k p p k+1 =: r k+1 p k+1 ,
where the degree of r k+1 is k+1 = k + -1. By induction,

d n dθ n r 1 p = r n+1 p n+1 ,
and r n+1 has degree n+1 = 1 + n( -1). Since m 11 is the quotient of the first order polynomial r 1 (s) = (2i -ξ 0 β -βs)/2, and the second order polynomial q(s), its n-th derivative, is

m (n) 11 (ξ 0 + θ) = r n+1 (ξ 0 + θ) q(ξ 0 + θ) n+1
and r n+1 is of degree n + 1. Using Lemma 5.1 and ( 5) we then obtain the required estimate,

|m (n) 11 (ξ 0 + θ)| ≤ C 1 + |ξ 0 + θ| n+1 (q 0 (1 + θ 2 )) n+1 ≤ C 1 + 2 n (|ξ 0 | n+1 + |θ| n+1 ) (q 0 (1 + θ 2 )) n+1 ≤ d n 1 + |θ| n+1 ,
where d n is independent of ξ 0 .

We are now ready to estimate the full phases φ a and φ g . Lemma 5.3. Let φ be either φ a or φ g . Suppose |δ| ≤ 1. Then there are constants c 0 and C n , independent of η, θ, ξ 0 and δ such that

|φ θ (δ, θ, η)| ≥ θ 2 16 , when |θ| ≥ c 0 1 + |η| 1/2 , ( 32
)
|∂ n θ φ(δ, θ, η)| ≤ C n      |θ| 2 + δ + |η|, n = 1, |θ| + δ, n = 2, 1 1+|θ| n-3 , n ≥ 3. ( 33
)
Additionally,

(34) φ(δ, θ, η) ≥ 0.

For φ a we have C n = 0 when n ≥ 4.

Proof. We first prove the statements for φ = φ a . Suppose |θ| ≥ c 0 (1 + |η|) 1/2 . Then |η| ≤ θ 2 /c 0 -1 and by ( 5)

|∂ θ φ a (δ, θ, η)| = |θ 2 + 2ηη 0 -δ| ≥ θ 2 -2|η| -1 ≥ θ 2 -2(θ 2 /c 2 0 -1) -1 ≥ (1 -2/c 2 0 )θ 2 , which gives (32) when c 0 ≥ 32/15. Similarly, |∂ θ φ a (δ, θ, η)| ≤ θ 2 + 2|η| + δ ≤ 2(θ 2 + |η| + δ),
showing (33) for n = 1. The bounds for larger n follow easily from an explicit calculation, yielding C n = 2 for 1 ≤ n ≤ 3 and C n = 0 for n ≥ 4.

To prove the statements for φ g we denote the correction term by w(δ, θ) = φ g -φ a . Lemma 5.2 gives

lim θ→±∞ w θ (δ, θ) θ 2 = lim θ→±∞ 1 θ 2 1 2 m 11 (ξ 0 + θ)(θ 2 -δ) 2 + 2m 11 (ξ 0 + θ)(θ 2 -δ)θ = - 1 4 + 1 = 3 4 .
Consequently, there is a K such that |w θ (δ, θ)| ≤ 7θ 2 /8 for all |θ| ≥ K, uniformly in ξ 0 and δ thanks to [START_REF] Zheng | Optimal error estimates for first-order Gaussian beam approximations to the Schrödinger equation[END_REF]. We now take c 0 = max(K, √ 32). Then for |θ| ≥ c 0 (1 + |η| 1/2 ) ≥ K, we have

|∂ θ φ g (δ, θ, η)| ≥ |∂ θ φ a (δ, θ, η)| -|∂ θ w(δ, θ)| ≥ (1 -2/c 2 0 )θ 2 - 7 8 θ 2 = (1/8 -2/c 2 0 )θ 2 ≥ 1 16 θ 2 ,
and (32) is proved.

For (33) we use (30) in Lemma 5.2. When n = 1 we have, as above,

|∂ θ φ g (δ, θ, η)| ≤ θ 2 + |δ -2η 0 η| + 1 2 |m 11 (ξ 0 + θ)|(θ 2 -δ) 2 + 2|m 11 (ξ 0 + θ)|θ(θ 2 -δ) ≤ θ 2 + |δ -2η 0 η| + d 1 2 |θ| 4 + δ 2 1 + θ 2 + 2d 0 |θ| 3 + δ|θ| 1 + |θ| ≤ 1 + 1 2 d 1 + 2d 0 (θ 2 + δ) + 2η 0 |η|,
which shows the result for n = 1 with C 1 = max(1 + d 1 /2 + 2d 0 , 2). For n = 2 we get

|∂ θθ φ g (δ, θ, η)| = -2θ + 1 2 m (2) 11 (ξ 0 + θ)(θ 2 -δ) 2 + 4m (1)
11

(ξ 0 + θ)(θ 2 -δ)θ + 2m 11 (ξ 0 + θ)(3θ 2 -δ) ≤ 2|θ| + d 2 2
|θ| 4 + δ 2 1 + |θ| 3 + 4d 1 |θ| 3 + δ|θ| 1 + θ 2 + 2d 0 3θ 2 + δ 1 + |θ| ≤ 2 + 1 2 d 2 + 4d 1 + 6d 0 (|θ| + δ) =: C 2 (|θ| + δ).
For n ≥ 3,

|∂ n θ w(δ, θ)| = 1 2 n =0 n m (n-) 11
(ξ 0 + θ) d dθ (θ 2 -δ) 2 = 1 2 max(4,n) =0 n m (n-) 11
(ξ 0 + θ) d dθ (θ 2 -δ) 2 ≤ C max(4,n) =0 n 1 1 + |θ| n-+1 (1 + |θ| 4-) ≤ C n 1 + |θ| n-3 , which shows the result for n ≥ 3 as ∂ 3 θ φ g = -1 + ∂ 3 θ w and ∂ n θ φ g = ∂ n θ w for n ≥ 4.
Finally, statement (34) for φ a is trivial, as φ a = 0, and for φ g it follows from (29) in Lemma 5.2 and (5), since

φ g (δ, θ, η) = 1 2 m 11 (ξ 0 + θ)(θ 2 -δ) 2 = (θ 2 -δ) 2 η 2 0 2|q(ξ 0 + θ)| 2 ≥ 0. This concludes the proof.
In the final part of this section we consider a space of function that is used in Lemma B.1. For a fixed phase function φ and order p we first introduce the basis functions

(35) W p (φ) = k φ (α k +1) φ : k α k = p, α k ≥ 1, ,
when p ≥ 1 and let W 0 (φ) be the constant function equal to one. Second, we denote by U p (φ) the linear span of these functions over the complex numbers, (36)

U p (φ) = span C W p (φ).
Functions in U p (φ) appear in Lemma B.1. Here we show that when φ is either φ a or φ g , these functions are bounded on subsets where the phase gradient grows at least quadratically.

Lemma 5.4.

Let K = {θ ∈ R | r 0 ≤ |θ| ≤ r 1 }, with 0 < R 0 ≤ r 0 < r 1 and, for some c > 0, |φ θ (δ, θ, η)| ≥ cθ 2 , for all θ ∈ K and |δ| ≤ 1,
where c and R 0 are independent of δ and η. Then for each u ∈ U p (φ(δ, •, η)), where φ is either φ a or φ g , there is a uniform bound |u(θ)| ≤ C, ∀θ ∈ K. The constant C depends on c and R 0 but is independent of r 0 , r 1 , δ and η.

Proof. We get from Lemma 5.3,

|φ θθ (δ, θ, η)| ≤ C 2 (|θ| + δ) ≤ C 2 (|θ| + 1) ≤ C 2 (1 + 1/R 0 )|θ|, ∀θ ∈ K,
while for n ≥ 3 we have

|∂ n θ φ(δ, θ, η)| ≤ C n 1 + |θ| n-3 ≤ C n , ∀θ ∈ R.
The constants C 2 and C n are independent of δ and η. Consider next w ∈ W p (φ(δ, •, η)) and assume it has M ≥ 1 factors and, without loss of generality, that the first M ≤ M factors have α k = 1. Then, for all θ ∈ K,

|w(θ)| = M k=1 |φ θθ (δ, θ, η)| |φ θ (δ, θ, η)| M k=M +1 |∂ α k +1 θ φ(δ, θ, η)| |φ θ (δ, θ, η)| ≤ M k=1 C 2 (1 + 1/R 0 )|θ| cθ 2 M k=M +1 C α k +1 |cθ 2 | ≤ max 2≤j≤p+1 C j M (1 + 1/R 0 ) M c M |θ| 2M -M ≤ C (1 + 1/R 0 ) M R 2M -M 0 =: C,
where C is independent of δ, η, r 0 and r 1 , but depends on c and R 0 . Since u ∈ U p (φ(δ, •, η)) is a linear combination of functions in W p (φ(δ, •, η)) the same bound holds true for u on K.

Estimates of oscillatory integrals

We consider integrals of the type

k 1/3 θ p r(θ)e ikφ(δ,θ,η) dθ,
where φ is either φ a or φ g and r ∈ W n,∞ (R), whose norm is defined by

||r|| W n,∞ (R) = n j=0 d j r dθ j L ∞ (R)
In general the integrand is then not in L 1 (R) and the integral must be defined in a generalized sense as an oscillatory integral. In this section, however, we only estimate the integral over bounded intervals that are defined using a smooth cutoff function ψ ∈ C ∞ c (R), which takes values in [0, 1], is equal to one on [-1, 1] and has supp(ψ) ⊂ [-2, 2]. This leaves us with integrals over compact domains with smooth integrands. Our main tool for estimating them are the identities stated in Lemma B.1. They are used to rewrite the integral on the domain where φ θ = 0. Lemma 5.3 in the previous section tells us when this is true. Lemma B.1 uses the space of functions U p defined in (36). Functions in U p are bounded on the domains we consider here, which is proved in Lemma 5.4. Together, Lemma B.1 and Lemma 5.4 constitute a precise version of the non-stationary phase lemma.

We will also make use of the simple inequalities

(37) ||uv|| W n,∞ (R) ≤ C n n j=0 j k=0 d k u dθ k L ∞ (R) d j-k v dθ j-k L ∞ (R) ≤ C n ||u|| W n,∞ (R) ||v|| W n,∞ (R) ,
for all u, v ∈ W n,∞ (R), and

(38) ||u(•/σ)|| W n,∞ (R) = n j=0 σ -j d j u dθ j L ∞ (R) ≤ max(1, σ -n )||u|| W n,∞ (R) , ∀u ∈ W n,∞ (R), σ = 0.
We start with an estimate of the integral between R and 2t where t is arbitrarily large. For this we consider a smooth cutoff around θ = R and θ = t > R and obtain bounds that are independent of t. Lemma 6.1. Let φ be either φ a or φ g and set

I t = k 1/3 (1 -ψ(θ/R))ψ(θ/t)θ p r(θ)e ikφ(δ,θ,η) dθ,
where c 0 ≤ R < t with c 0 as in Lemma 5.3, |δ| ≤ 1 and r ∈ W n,∞ (R). If c 0 (1 + |η| 1/2 ) ≤ R and n ≥ 1 + p/2, there is a constant C n independent of k, R, δ, η and t such that

|I t | ≤ C n k 1/3-n ||r|| W n,∞ (R) .
Proof. On this domain the results in Section 5 show that the phase gradient does not vanish, and |φ θ | ≥ cθ 2 . Since the integrand is smooth and compactly supported we can therefore use the non-stationary phase lemma to estimate the integral. For sufficiently regular r the repeated partial integrations in this lemma enables us to offset the growing θ p factor and obtain a bound that is independent of t.

To be precise, let

b(θ) = (1 -ψ(θ/R))ψ(θ/t)r(θ), which is supported in the compact set K = {θ ∈ R | R ≤ |θ| ≤ 2t}.
Then by Lemma 5.3 we have |φ θ (δ, θ, η)| ≥ θ 2 /16 on K, independent of δ, η and t, since |θ| ≥ R > c 0 (1 + |η| 1/2 ). We apply Lemma B.1 with a(θ) = b(θ)θ p and let D be any bounded open set containing K. This gives

I t = k 1/3 D b(θ)θ p e ikφ(δ,θ,η) dθ = k 1/3 (ik) -n n =0 K d dθ b(θ)θ p u ,n (θ) φ θ (δ, θ, η) n e ikφ(δ,θ,η) dθ, u ,n ∈ U n-(φ(δ, •, η)). ( 39
)
where the space U p is defined in (36). Since K satisifies all conditions in Lemma 5.4 and R ≥ c 0 > 0 we obtain a uniform bound,

|u ,n (θ)| ≤ C ,n , ∀θ ∈ K,
where C ,n depends on c 0 , but is independent of δ, η, R and t. This allows us to estimate I t as

|I t | ≤ Ck 1/3-n n =0 min( ,p) j=0 K |b ( -j) (θ)| |θ| p-j |θ| 2n dθ ≤ Ck 1/3-n n =0 min( ,p) j=0 ||b ( -j) || L ∞ (R) ∞ R dθ |θ| 2n-p+j ≤ Ck 1/3-n n =0 ||b ( ) || L ∞ (R) p 1 R 2n-p+j+1 ≤ C k 1/3-n c 2n-p+1 0 ||b|| W n,∞ (R) ,
where we also used the fact that 2(n -1) ≥ p and c 0 ≤ R ≤ t. Moreover, by (37) and (38),

||b|| W n,∞ (R) ≤ C 2 n max(1, c -n 0 ) 2 ||1 -ψ|| W n,∞ (R) ||ψ|| W n,∞ (R) ||r|| W n,∞ (R) ≤ C||r|| W n,∞ (R)
The result in the lemma follows.

Next we consider the main part of the integral for small η and δ with the Airy phase. The estimate involves the norm of r with an argument scaled by k 1/3 . Lemma 6.2. Let I = k 1/3 ψ(θ/R)θ p r(θ)e ikφa(δ,θ,η) dθ, where 0 < R 0 ≤ R ≤ R 1 and r ∈ W n,∞ (R) with n = (p + 1)/5 . Moreover, suppose

(40) 2|η| + |δ| ≤ ξ 2 0 2 k -2/3 , k ≥ 1.
Then there is a constant C, depending on R 0 and R 1 , but independent of k, R, η and δ such that

|I| ≤ Ck -p/3 ||r k || W n,∞ (R) , rk (ζ) := r k -1/3 ζ ,
Proof. For δ -2ηη 0 = 0 the phase φ a has a degenerate stationary point at the origin. We will therefore treat the integral in the vicinity of the origin separately. Away from the origin we have the same type of lower bound ∂ θ φ a ≥ cθ 2 as in Lemma 6.1 and we can therefore once again use the non-stationary phase lemma to estimate the integral.

For the proof we use the rescaled variables

ζ = k 1/3 θ, δ = k 2/3 δ, η = k 2/3 η. Since kφ a ( δ/k 2/3 , ζ/k 1/3 , η/k 2/3 ) = - 1 3 ζ 3 + ζ( δ -2η 0 η) = φ a ( δ, ζ, η),
we can rewrite the integral as

I = k -p/3 ψ ζ k 1/3 R ζ p rk (ζ)e iφa( δ,ζ,η) dζ.
We then divide the integral into two pieces,

I = k -p/3 ψ (ζ) ψ ζ k 1/3 R ζ p rk (ζ)e iφa( δ,ζ,η) dζ + k -p/3 (1 -ψ (ζ))ψ ζ k 1/3 R ζ p rk (ζ)e iφa( δ,ζ,η) dζ = k -p/3 (I 1 + I 2 ),
where I 1 is the part close to the origin containing the stationary point, and I 2 is the remaining part. Note that I 2 matches the general form of the integral in Lemma 6.1 if we take t = k 1/3 R 2 .

For I 1 we simply have

|I 1 | ≤ 2 -2 |ζ| p |r k (ζ)| dζ ≤ C||r k || L ∞ (R) ,
with C independent of k. If 2k 1/3 R ≤ 1 we have I 2 = 0 and the proof is complete. We assume henceforth that 2k 1/3 R > 1 and let

b(ζ) = (1 -ψ(ζ))ψ ζ k 1/3 R rk (ζ),
the support of which lies in the compact set

K = {ζ ∈ R | 1 ≤ |ζ| ≤ 2k 1/3 R}. Then for ζ ∈ K, by (40), |∂ ζ φ a ( δ, ζ, η)| = |ζ 2 + k 2/3 (2η 0 η -δ)| ≥ |ζ 2 | -k 2/3 (2|η| + |δ|) ≥ |ζ 2 | - ξ 2 0 2 ≥ 1 2 |ζ 2 |.
I 2 = (ik) -n n =0 K d dζ b(ζ)ζ p u ,n (ζ) φζ ( δ, ζ, η) n e ikφa( δ,ζ,η) dζ,
where u ,n ∈ U n-, with U p defined in (36). This expression is now estimated in the same way as (39) above. Since | δ| ≤ ξ 2 0 /2 ≤ 1/2 and K satsfies the assumptions of Lemma 5.4 we obtain a uniform bound for u ,n on K. Then

|I 2 | ≤ Ck -n n =0 min( ,p) j=0 ||b ( -j) || L ∞ (R) 2k 1/3 R 1 dζ |ζ| 2n-p+j ≤ Ck -n ||b|| W n,∞ (R) 2k 1/3 R 1 dζ |ζ| 2n-p ≤ Ck -n ||b|| W n,∞ (R) max(1, (2k 1/3 R) p-2n )2k 1/3 R ≤ C||b|| W n,∞ (R) max(k 1/3-n , k (p-5n+1)/3 R p-2n )R 1 ≤ C max(1, R p-2n 0 , R p-2n 1 )R 1 ||b|| W n,∞ (R) =: C||b|| W n,∞ (R) ,
since p + 1 ≤ 5n and k ≥ 1. Moreover, as in the proof of Lemma 6.1, by ( 37) and (38), since

k 1/3 R > 1/2, |I 2 | ≤ C||b|| W n,∞ (R) ≤ CC 2 n 2 n ||1 -ψ|| W n,∞ (R) ||ψ|| W n,∞ (R) ||r k || W n,∞ (R) ≤ C||r k || W n,∞ (R)
. Together the estimates of I 1 and I 2 then prove the lemma. We finally note that since r ∈ W n,∞ (R) the norm ||r k || W n,∞ (R) is bounded because of (38) and (40).

Finally, we show that the derivatives of the Airy function are well approximated by an oscillatory integral with a monomial factor and the Airy phase.

I - 2πi p k p/3 Ai (p) (k 2/3 (2η 0 η -δ)) ≤ C n k 1/3-n .
Proof. The Fourier transform of Ai (p) in S (R) is

F(Ai (p) )(ζ) = 1 √ 2π (iζ) p e i ζ 3 3 .
Therefore, using Lemma 2.1, and noting that

F -1 (•)(ρ) = F(•)(-ρ), Ai (p) (ρ) = F -1 1 √ 2π (iζ) p e i ζ 3 3 (ρ) = lim t→∞ 1 2π ψ(ζ/t)(iζ) p e i ζ 3
3 +ζρ dζ.

After rescaling θ = k -1/3 ζ we get

ψ(ζ/t)(iζ) p e i ζ 3 3 +ζρ dζ = i p k p+1 3 ψ(θk 1/3 /t)θ p e ik θ 3
3 +k -2/3 θρ dθ.

It follows that if

(41) ρ = k 2/3 (2η 0 η -δ), then 2π (ik 1/3 ) p Ai (p) (k 2/3 (2η 0 η -δ)) = k 1/3 lim t→∞ ψ(θk 1/3 /t)θ p e ik θ 3 3 +k -2/3 θρ dθ = k 1/3 lim t→∞ ψ(θ/t)θ p e -ikφa(δ,θ,η) dθ = k 1/3 (-1) p lim t→∞ ψ(-θ/t)θ p e -ikφa(δ,-θ,η) dθ.
Since φ a is odd and ψ is even in θ we obtain

2πi p k p/3 Ai (p) (k 2/3 (2η 0 η -δ)) = I + k 1/3 lim t→∞ (1 -ψ(θ/R))ψ(θ/t)θ p e ikφa(δ,θ,η) dθ,
The result now follows from Lemma 6.1, with r ≡ 1.

Solution estimates

In Sections 3 and 4 it was shown that the partial Fourier transform in y for both the exact solution and the Gaussian beam approximation can be written on the form

û(x, η) = v(x, η; k) Â(η),
where A is the amplitude function and v for the two cases are given in [START_REF] Lafitte | The kernel of the Neumann operator for a strictly diffractive analytic problem[END_REF] and (25). In this section we prove bounds of those v in terms of k and η, which are valid for all η ∈ R, x ∈ [0, x c ] and k ≥ 1. We start with the Gaussian beam superposition case and estimate vGB as follows.

Lemma 7.1. For vGB defined in (25), there is a constant M such that

(42) |v GB (η, x, k)| ≤ M 1 + log 1 + |η| 1/2 k 1/2 , for all η ∈ R, x ∈ [0, x c ] and k ≥ 1.
Proof. From (24) in Section 4.2 we have

I t = k 1/3 ψ(θ/t)r(θ)e ikφg(δ,θ,η) dθ, r(θ) = 1 2π
q(ξ 0 ) q(ξ 0 + θ) .

We note that r ∈ W n,∞ (R) for all n by Lemma 5.1 and that |δ|

= |x c -x| ≤ |x c | = ξ 2 0 ≤ 1.
We divide the integral into two parts. Let R = R(η) = c 0 1 + |η| 1/2 with c 0 as in Lemma 5.3 and define

I t = k 1/3 ψ(θ/R)ψ(θ/t)r(θ)e ikφg(δ,θ,η) dθ + k 1/3 (1 -ψ(θ/R))ψ(θ/t)r(θ)e ikφg(δ,θ,η) dθ =: I 1,t + I 2,t .
For I 1,t we have, again using Lemma 5.1, and that fact that φ g ≥ 0 by Lemma 5.3,

|I 1,t | ≤ k 1/3 q 1 q 0 2R(η) -2R(η) e -k φg(x,θ,η) √ 1 + θ 2 dθ ≤ k 1/3 q 1 q 0 2R(η) -2R(η) dθ √ 1 + θ 2 ≤ Ck 1/3 log(R(η)) ≤ C k 1/3 log(1 + |η| 1/2 ).
For I 2,t we use Lemma 6.1 with p = 0, which says that for any n ≥ 1 and t > R(η),

|I 2,t | ≤ C n k 1/3-n ||r|| W n,∞ (R) ≤ C n k 1/3-n ,
where C n is independent of k, R(η), δ, η and t. By (23) we have

|P GB | = 2(πξ 0 ) 1/2 ≤ 2 √
π for all η and k. Then we get

|v GB (η, x, k)| ≤ lim t→∞ k 1/6 |P GB (x, η)|(|I 1,t | + |I 2,t |) ≤ 2 √ πk 1/6 C k 1/3 log(1 + |η| 1/2 ) + C n k 1/3-n ,
and upon taking n = 1 the Lemma follows with M = 2 √ π max(C , C 1 ).

Next, for the exact solution, we estimate v.

Lemma 7.2. For v defined in [START_REF] Lafitte | The kernel of the Neumann operator for a strictly diffractive analytic problem[END_REF] there is a constant M such that

|v(η, x, k)| ≤ M k 1/6 , for all η ∈ R, x ∈ [0, x c ] and k ≥ 1.
Proof. From [START_REF] Lafitte | The kernel of the Neumann operator for a strictly diffractive analytic problem[END_REF] we have

v(η, x, k) = k 1/6 P (k, η)Ai(k 2 3 (x -X)), P (k, η) = ᾱk -1/6 Ai(αk 2/3 X) , X(η) = 1 -(η 0 + η) 2 .
Then, using ( 68) and (70) in Lemma 9.1,

|v(η, x, k)| = |Ai(k 2 3 (x -X)| |Ai(αk 2 3 X)| ≤ C (1 + k 2 3 |x -X|) -1/4 (1 + k 2 3 |X|) -1/4 ≤ C 1 + k 2 3 (|x| + |x -X|) 1 + k 2 3 |x -X| 1/4 = C 1 + k 2 3 |x| 1 + k 2 3 |x -X| 1/4 ≤ C 1 + |x c |k 2 3 1/4 ≤ C(1 + k 1/6 ).
The result follows as k ≥ 1. (The estimate is sharp for x = X.)

Proof of the Main Result

In this section we prove the main result Theorem 1.1 estimating the L ∞ error between the exact solution and the Gaussian beam solution. To estimate the difference between u GB (x, y) and the exact solution u(x, y), it is enough to control the L 1 norm of the difference between their scaled Fourier transforms since

(43) ||u GB (x, •) -u(x, •)|| L ∞ ≤ k 2π ||û GB (x, •) -û(x, •)|| L 1 .
We will use this strategy. From ( 25) and ( 10) we get ûGB (x,

η 0 + η) -û(x, η 0 + η) = vGB (η, x, k) -v(η, x, k) Â(η).
We divide the expression into two parts, one for |η| smaller than O(k -2/3 ) and one for |η| larger than O(k -2/3 ). Thus, for c to be determined below, we let

k 2π ||û GB (x, •) -û(x, •)|| L 1 ≤ k 2π |η|≤ck -2/3 |v GB (x, η) -v(x, η)|| Â(η)|dη (44) + k 2π |η|≥ck -2/3 |v GB (x, η) -v(x, η)|| Â(η)|dη =: E 1 + E 2 .
For the large values of |η| we can immediately get a bound of O(k -1 ) by using the fact that  has very rapid decay, being the Fourier transform of A ∈ S. This is used in the following lemma.

Lemma 8.1. Suppose f (•, k) ∈ L 1 loc (R) for each k ≥ 1 and A ∈ S(R). Let c, β > 0 be given. If there exist r ∈ R and q, M ∈ R + such that,

(45) |f (η, k)| ≤ M (1 + |kη| q )k r , ∀k ≥ 1, |kη| ≥ ck β ,
then, for each p ≥ 0, there exists a constant C p , independent of k, but dependent on A, M, q, r, c, such that

k 1/2 |kη|≥ck β f (η, k) Â(η) dη ≤ C p k -p .
Proof. Since A ∈ S(R) for all ≥ 0 there exists c such that

|F(A)(η)| ≤ c /(1 + |η| ). Hence, | Â(η)| = | √ kF(A)(kη)| ≤ c √ k 1 + |kη| , ∀η.
Then for > q + 1, with ξ = kη,

k 1/2 |kη|≥ck β f (η, k) Â(η) dη ≤ c k |kη|≥ck β |f (η, k)| 1 + |kη| dη ≤ c M k r+1 |kη|≥ck β 1 + |kη| q 1 + |kη| dη = c M k r |ξ|≥ck β 1 + |ξ| q 1 + |ξ| dξ ≤ c M k r |ξ|≥ck β 1 + |ξ| q |ξ| dξ = 2M c c 1- -1 k r+β(1-) + c 1-+q -1 -q k r+β(1-+q) .
For the given p we now take = p := max((r + p)/β + q + 1, q + 2). Then r + β(1 -+ q) ≤ -p and

k 1/2 |kη|≥ck β f (η, k) Â(η) dη ≤ 2M c p c 1-p p -1 k -p-βq + c 1-p +q p -1 -q k -p ≤ 2M c p c 1-p p -1 + c 1-p +q p -1 -q =:Cp k -p ,
which is the desired estimate.

In our case we let f = (v GB -v)/ √ 2π for fixed x. Then it follows from Lemma 7.1 and Lemma 7.2 that f satisfies (45), as for k ≥ 1,

|v GB (η, x, k)-v(η, x, k)| ≤ M k 1/6 +M (1+log(1+|η| 1/2 ))k 1/2 ≤ M (1+|η| 1/2 )k 1/2 ≤ M (1+|kη| 1/2 )k 1/2 ,
where M is in fact independent also of x ∈ [0, x c ]. Lemma 8.1 with β = 1/3, q = 1/2, r = 1/2, c = ξ 2 0 /4 and p = 1 now shows that

(46) |E 2 | ≤ Ck -1 .
The choice of c will be motivated below in the next step.

To estimate E 1 we will use more precise estimates of |v GB (η, x, k) -v(η, x, k)| for small |η| and the following lemma.

Lemma 8.2. Suppose f (•, k) ∈ L 1 loc (R) for each k ≥ 1 and A ∈ S(R). Let c, β > 0 be given. If there exist r ∈ R and q, M ∈ R + such that, (47) |f (η, k)| ≤ M (1 + |kη| q )k r , ∀k ≥ 1, |kη| ≤ ck β ,
then there exists a constant C, independent of k, but dependent on A, M, q, r, c, such that

k 1/2 |kη|≤ck β f (η, k) Â(η) dη ≤ Ck r .
Proof. As in the proof of Lemma 8.1, when > q + 1 we get

k 1/2 |kη|≤ck β f (η, k) Â(η) dη ≤ c M k r |ξ|≤ck β 1 + |ξ| q 1 + |ξ| dξ ≤ c M 1 + |ξ| q 1 + |ξ| dξ =:C k r .
This proves the lemma.

As for E 2 we apply this lemma with f = (v GB -v)/ √ 2π and to get the bound (47) we need to estimate the difference between vGB and v when |kη| ≤ ck 1/3 . This estimate is the main part of the proof.

To examine vGB -v more carefully we first recall the expressions:

vGB (η, x, k) = k 1/6 P GB (k, η)I(η, x, k), v(η, x, k) = k 1/6 P (k, η)Ai(k 2 3 (x -X(η))),
where

I(η, x, k) = lim t→∞ k 1/3 q(ξ 0 ) 2π ψ(θ/t) e ikφg(x,θ,η) q(ξ 0 + θ) dθ P GB (k, η) = 2(πξ 0 ) 1/2 e -iπ/4 e 2 3 ikξ 3 0 -2ikη0ξ0η , and (48)
X(η) := 1 -(η + η 0 ) 2 , P (k, η) = ᾱk -1/6 Ai(αk 2/3 X) .
We divide the difference vGB (η,

x, k) -v(η, x, k) into three parts v(η, x, k) -vGB (η, x, k) = k 1/6 P GB (k, η)[Ai(k 2 3 (x -X)) -Ai(k 2 3 (x -X -η 2 ))] + k 1/6 [P (k, η) -P GB (k, η)]Ai(k 2 3 (x -X)) + k 1/6 P GB (k, η)[Ai(k 2 3 (x -X -η 2 )) -I(η, x, k)] =: R 1 + R 2 + R 3 .
In three Lemmas below we show that when x = x c and η ≤ ξ 2 0 k -2/3 /4, there is a constant M such that

|R 1 | ≤ M |kη| 2 k -1 , |R 2 | ≤ M (1 + |kη| 2 )k -5/6 , |R 3 | ≤ M (1 + |kη| 2 )k -5/6 .
It follows that

|v(η, x c , k) -vGB (η, x c , k)| ≤ |R 1 + R 2 + R 3 | ≤ 3M (1 + |kη| 2 )k -5/6 , when |kη| ≤ ξ 2 0 4 k 1/3 and k ≥ 1,
Then applying Lemma 8.2 with β = 1/3, q = 2, r = -5/6 and c = ξ 2 0 /4 shows that |E 1 | ≤ Ck -5/6 , when x = x c . Together with (43), ( 44) and (46) this proves Theorem 1.1.

Note that the estimates of R 1 and R 2 above are shown to be valid for all x ∈ [0, x c ], while the R 3 estimate is considered, in this paper, only for x = x c . Furtheremore, note that R 2 and R 3 exhibit the same loss of decay through the term k -5/6 . In R 2 this comes from the estimate (52) and R 3 has k 1/6 embedded in (54).

We now turn to proving the lemmas about R j .

8.1. Estimate of R 1 .

Lemma 8.3. There is a constant M independent of η and k ≥ 1, such that

|R 1 | ≤ M kη 2 , when |η| ≤ 1. Proof. Since |P GB | ≤ 2(πξ 0 ) 1/2 we have |R 1 | ≤ 2 √ πk 1/6 |Ai(k 2 3 (x -X)) -Ai(k 2 3 (x -X -η 2 ))|.
Moreover, from (69) in Lemma 9.1 we get

|Ai(k 2 3 (x -X) -Ai(k 2 3 (x -X -η 2 )| = k 2 3 η 2 0 Ai (k 2 3 (x -X -s))ds ≤ C 3 k 2 3 η 2 max 0≤s≤η 2 (1 + |k 2 3 (x -X -s)|) 1/4 . Then, since |X| ≤ (|η| + |η 0 |) 2 -1 ≤ 3, |x| ≤ 1, s ≤ η 2 ≤ 1, we obtain |Ai(k 2 3 (x -X)) -Ai(k 2 3 (x -X -η 2 ))| ≤ C 3 k 2 3 η 2 (1 + 5k 2
|R 2 | ≤ M (1 + k 2 η 2 )k r , when |η| ≤ ξ 2 0 4 k -2/3 , r = -5/6, x = x c , -1, 0 ≤ x < x c .
Proof. We start by proving two estimates of X(η). We use the inequalities 1 -x ≤ √ 1 -x ≤ 1 -x/2 which hold for x ∈ [0, 1]. The definition (48) together with the assumption on η and the fact that k ≥ 1, then gives

X(η) = 1 -(η + η 0 ) 2 ≥ 1 - ξ 2 0 4 + 1 -ξ 2 0 2 ≥ 1 -1 - ξ 2 0 4 2 ≥ ξ 2 0 4 =: X 0 > 0.
Moreover, since

x c = ξ 2 0 , |X(η) -x c | = |η||2η 0 -η| ≤ |η| 2 1 -ξ 2 0 + ξ 2 0 4 ≤ |η| 2 -ξ 2 0 + ξ 2 0 4 ≤ 2|η|.
Clearly, we also have X(η) ≤ 1, and therefore, in summary,

(49) 0 < X 0 ≤ X(η) ≤ 1, |X(η) -x c | ≤ 2|η|.
Next, we rewrite P GB , adopting the the defintion and

(50) Ai(z) = 1 2 √ π z -1 4 e -
P GB = ᾱk -1/6
Ai(αk 2/3 x c )

e -2ikηη0 √ xc = ᾱk -1/6
Ai(αk 2/3 X)

x c X 1/4 e -2 3 ik(X 3/2 -x 3/2 c )-2ikηη0 √ xc =: ᾱk -1/6 g(η)
Ai(αk 2/3 X) .

We get

k 1/6 |P -P GB | = 1 Ai(αk 2/3 X) - g(η)
Ai(αk 2/3 X)

≤ 1 Ai(αk 2/3 X) - 1
Ai(αk 2/3 X)

+ 1 -g(η)
Ai(αk 2/3 X)

= 1 | Ai(αk 2/3 X)| Ai(αk 2/3 X) -Ai(αk 2/3 X) Ai(αk 2/3 X) + |1 -g(η)| .
We can then estimate R 2 as

|R 2 | ≤ |Ai(k 2 3 (x -X))| | Ai(αk 2/3 X)| Ai(αk 2/3 X) -Ai(αk 2/3 X) Ai(αk 2/3 X) + |1 -g(η)| . ( 51
)
We will now study the different parts of this expression separately.

• Estimate of Ai(αk 2/3 X)-Ai(αk 2/3 X) Ai(αk 2/3 X) . This is given directly by (67) in Lemma 9.1 with s 0 = X 0 , as then k 2/3 X(η) ≥ k 2/3 X 0 ≥ s 0 . We get Ai(αk 2/3 X) -Ai(αk 2/3 X)

Ai(αk 2/3 X)

≤ C 1 |k 2/3 X| -3/2 ≤ C 1 X -3/2 0 k -1 =: D 1 k -1 ,
where C 1 is the constant in (67).

• Estimate of |1 -g(η)|. Using Taylor's formula for x → x 3/2 around x = X(0) = 1 -η 2 0 = ξ 2 0 = x c , we compute

X(η) 3/2 = X(0) 3/2 + 3 2 X(0) 1/2 (X(η) -X(0)) + R(X(η) -X(0)) 2 , = x 3/2 c + 3 2 √ x c (2ηη 0 -η 2 ) + R(X(η) -X(0)) 2 , |R| ≤ sup ξ≥X0 3 8 ξ -1/2 = 3 8 √ X 0 . Therefore, 2 3 (X 3/2 -x 3/2 c ) + 2ηη 0 √ x c = - √ x c η 2 + 2 3 R(X(η) -x c ) 2 ,
and consequently, by (49),

2 3 (X 3/2 -x 3/2 c ) + 2ηη 0 √ x c ≤ |η| 2 + |X(η) -x c | 2 4 √ X 0 ≤ |η| 2 + |η| 2 √ X 0 =: D 2 η 2 .
This gives us

|1 -g(η)| ≤ 1 -e -ik( 2 3 (X 3/2 -x 3/2 c )+2ηη0 √ xc) + x c X(η) 1/4 -1 ≤ k 2 3 (X 3/2 -x 3/2 c ) + 2ηη 0 √ x c + 1 + 2η X 0 1/4 -1 ≤ D 2 kη 2 + 1 2X 0 η ≤ D 3 (1 + k 2 η 2 )k -1 , for D 3 = max(D 2 , (2X 0 ) -1 ). • Estimate of |Ai(k 2 3 (x-X))| | Ai(αk 2/3 X)| .
We divide this into three subcases. Suppose first that

|x -x c | > |X(η) -x c |.
Then by (49),

k 2/3 |x -X(η)| ≥ k 2/3 |x -x c | -k 2/3 |X(η) -x c | ≥ k 2/3 |x -x c | -2k 2/3 |η| ≥ k 2/3 |x -x c | - ξ 2 0 2 .
By (68) in Lemma 9.1 and (50),

|Ai(k 2 3 (x -X))| | Ai(αk 2/3 X)| ≤ 2 √ πk 1/6 |X| 1/4 C 2 (1 + |k 2 3 (x -X)|) -1/4 ≤ 2 √ πk 1/6 C 2 (1 -ξ 2 0 /2 + k 2 3 |x -x c |) -1/4 ≤ 2 √ πC 2 |x -x c | -1/4 ,
where C 2 is the constant in (68). On the other hand, if

|X(η) -x c | ≥ |x -x c | > 0,
then by (49),

k 2/3 ≤ ξ 2 0 4|η| ≤ ξ 2 0 2|X(η) -x c | ≤ ξ 2 0 2|x -x c | ,
and we obtain the same estimate as above, via

|Ai(k 2 3 (x -X))| | Ai(αk 2/3 X)| ≤ 2 √ πk 1/6 |X| 1/4 C 2 (1 + |k 2 3 (x -X)|) -1/4 ≤ 2 √ πk 1/6 C 2 ≤ 2 √ πC 2 ξ 2 0 2|x -x c | 1/4 ≤ 2 √ πC 2 |x -x c | -1/4 .
Finally, when x = x c (caustic case) we can not get better than

|Ai(k 2 3 (x -X))| | Ai(αk 2/3 X)| ≤ 2 √ πk 1/6 |X| 1/4 C 2 (1 + |k 2 3 (x -X)|) -1/4 ≤ 2 √ πC 2 k 1/6 .
In summary, we have with

D 4 = 2 √ πC 2 , (52) |Ai(k 2 3 (x -X))| | Ai(αk 2/3 X)| ≤ D 4 k 1/6 , x = x c , |x -x c | -1/4 , x < x c .
We can now put the estimates together and apply them to R 2 in (51). We get

|R 2 | ≤ D 1 k -1 + D 3 (1 + k 2 η 2 )k -1 D 4 k 1/6 , x = x c , |x -x c | -1/4 , x < x c , ≤ M (1 + k 2 η 2 )k r ,
where

M = max(D 1 , D 3 )D 4 1, x = x c , |x -x c | -1/4 , x < x c .
This proves the lemma. 8.3. Estimate of R 3 at the caustic. This is the main estimate. Here we assume that x = x c . Lemma 8.5. For x = x c there is a constant M independent of η and k ≥ 1, such that

|R 3 | ≤ M (1 + k 2 η 2 )k -5/6 , when |η| ≤ ξ 2 0 4 k -2/3 .
Proof. We consider ρ = 2k 2/3 η 0 η, which amounts to taking δ = 0 in (41). By the assumption on η and k it is bounded as

(53) |ρ| ≤ 2k 2/3 |η| ≤ ξ 2 0 /2 ≤ 1 2 .
Moreover, since k 2/3 (x c -X -η 2 ) = ρ and as before, |P GB | ≤ 2(πξ 0 ) 1/2 , we get (54)

|R 3 | = k 1/6 |P GB (k, η)| • |Ai(k 2 3 (x c -X -η 2 )) -I(η, x c , k)| ≤ Ck 1/6 |Ai(ρ) -I(η, x c , k)|.
Hence, we need to estimate |Ai(ρ) -I|. Let (55) r(θ) = 1 2π q(ξ 0 ) q(ξ 0 + θ) .

As in the proof of Lemma 7.1 we then use the fact that I = lim t→∞ I t where I t is defined and divided as

I t = k 1/3 ψ(θ/t)r(θ)e ikφg(0,θ,η) dθ = k 1/3 ψ(θ/R)ψ(θ/t)r(θ)e ikφg(0,θ,η) dθ + k 1/3 (1 -ψ(θ/R))ψ(θ/t)r(θ)e ikφg(0,θ,η) dθ =: I main + I t,tail .
With c 0 as in Lemma 5.3 we choose here R = 3c 0 /2, independent of η, which implies that for all η which we consider, (56)

c 0 (1 + |η| 1/2 ) ≤ c 0 (1 + ξ 0 /2) ≤ 3 2 c 0 = R.
Moreover, we take t ≥ 2R = 3c 0 such that ψ(θ/R)ψ(θ/t) = ψ(θ/R). To analyze I main we then first note that it can be written as

I main = k 1/3 ψ(θ/R)r(θ)e ikφa(0,θ,η) dθ, r(θ) = r(θ)e ik 1 2 m11(ξ0+θ)θ 4 .
We next expand r in terms of θ, first using the Taylor expansion of exp(iz),

e iz = 1 + iz + (iz) 2 2 + (iz) 3 Z(z), Z(z) = 1 2i 1 0 e isz (1 -s) 2 ds. This gives r(θ) = r(θ) 1 + ik 2 m 11 (ξ 0 + θ)θ 4 + (ik) 2 8 m 11 (ξ 0 + θ) 2 θ 8 + (ik) 3 8 r(θ)Z k 2 m 11 (ξ 0 + θ)θ 4 m 11 (ξ 0 + θ) 3 θ 12 .
Furthermore, let v (θ) := r(θ)m 11 (ξ 0 + θ) ,

and Taylor expand these functions as

v (θ) = p j=0 v (j) (0)θ j j! + V p, (θ)θ p+1 , V p, (θ) = 1 p! 1 0 v (p+1) (tθ)(1 -t) p dt. Then r(θ) = v 0 (0) + v 0 θ + 1 2 v 0 (0)θ 2 + V 2,0 (θ)θ 3 + ik 2 v 1 (0)θ 4 + ik 2 v 1 (0)θ 5 + ik 2 V 1,1 (θ)θ 6 + (ik) 2 8 v 2 (0)θ 8 + (ik) 2 8 V 0,2 (θ)θ 9 + (ik) 3 8 v 3 (θ)Z k 1 2 m 11 (ξ 0 + θ)θ 4 θ 12 .
From this expansion of r we now get a corresponding expansion of I main ,

I main = I S + I V2,0 + 1 2 I V1,1 + 1 8 I V0,2 + I Z , (57)
I S = v 0 (0)I 0 + v 0 (0)I 1 + 1 2 v 0 (0)I 2 + v 1 (0) ik 2 I 4 + ik 2 v 1 (0)I 5 + (ik) 2 8 v 2 (0)I 8 ,
where I p = k 1/3 ψ(θ/R)θ p e ikφa(0,θ,η) dθ, I V p, = (ik) k 1/3 ψ(θ/R)V p, (θ)θ 3 +3 e ikφa(0,θ,η) dθ, and

I Z = (ik) 3 8 k 1/3 ψ(θ/R)w(θ, k)θ 12 e ikφa(0,θ,η) dθ, with w(θ, k) = v 3 (θ)z(θ, k), z(θ, k) = Z k 2 m 11 (ξ 0 + θ)θ 4 .
We will next show that the last four terms in (57) are at most of size O(1/k). To see this, we note that by Lemma 5.1 and Lemma 5.2, both r and m 11 (ξ 0 + •) belong to W n,∞ (R) for all n, ξ 0 , and their W n,∞ -norms are bounded independent of ξ 0 . By (37) the same is true for v , for all . Therefore, by (38),

||V p, || W n,∞ (R) ≤ 1 p! 1 0 ||v (p+1) (• t)|| W n,∞ (R) (1 -t) p dt ≤ 1 p! 1 0 max(1, t n )||v (p+1) || W n,∞ (R) (1 -t) p dt ≤ 1 p! ||v || W n+p+1,∞ (R) ,
showing that also V p, ∈ W n,∞ (R) for all n, p, . Since (40) is satisfied under the assumptions on η, δ and k, we can use Lemma 6.2 with n = (3 + 4)/5 to estimate

|I V p, | ≤ Ck -(3 +3)/3 ||V p, (•/k 1 3 )|| W n,∞ (R) ≤ Ck -1 max(1, k -n 3 )||V p, || W n,∞ (R) ≤ Ck -1 .
For I Z we first observe that

z(θ/k 1 3 , k) = Z k -1 3 2 m 11 (ξ 0 + k -1 3 θ)θ 4 .
It remains to estimate I t,tail . By (56) we get from Lemma 6.1 with p = 0 and n = 2, for all t > R, that (61)

|I t,tail | ≤ C 2 k 1/3-2 ||r|| W 2,∞ (R) ≤ Ck -5/3 ,
where the constant is independent of t. In conclusion, using (58), ( 59), ( 60) and (61) we have shown that

|I(η, x c , k) -Ai(ρ)| ≤ |I main -I S | + |I S -ĨS | + | ĨS -Ai(ρ)| + lim t→∞ |I t,tail | ≤ C 1 + |kη| 2 + |kη| + k -2/3 k -1 ≤ C 1 + |kη| 2 k -1 .
Together with (54) this concludes the proof of Lemma 8.5. Finally note that, away from the caustic point, i.e. x < x c , the method used here to estimate R 3 will not give sharp results; if the stationary phase method is applied directly to φ g extra decay in k follows.

Properties of the Airy function

Here we show some known properties of the Airy function and we derive a few consequences in two lemmas. A more complete source for information about Airy functions is [12], which we frequently cite below. We consider the Airy function of the first kind Ai and second kind Bi.

(P1) The Airy functions are linearly independent solutions of the Airy differential equation Then, for real s > 0,

Ai(-s) = 1 2 π -1 2 s -1 4 cos 2 3 s 3 2 - π 4 + i sin 2 3 s 3 2 - π 4 ,
and it follows easily from [12, Section 9.7 (ii,iii)] that

(64) Ai(s) -Ai(s) ≤ Cs -3/2 Ai(s) , Ai(-s) -2 Ai(-s) ≤ Cs -3/2 (65) Ai (s) + √ s Ai(s) ≤ Cs -1 Ai(s) , Ai (-s) -2 √ s Ai(-s) ≤ Cs -1 .
We can now prove the following lemmas. Ai(βs) s 3/2 . Now suppose s ≥ s 0 > 0. We first note that for all s > 0,

|Φ 0 (βs)| ≥ 1 -|Φ 0 (βs) -1| ≥ 1 - |c 1 | s 3/2 .
Hence, for s ≥ By (P2) we also have |Ai(αs)| ≥ c for 0 ≤ s ≤ s 0 and some c > 0. Moreover, the identity [12, Eq. 9. where p m and q m are polynomials given by the recursions (73) p m+1 = p m + xq m , q m+1 = p m + q m , p 0 = 1, q 0 = 0. Furthermore, (74) Ai (3p+2) (0) = 0, p = 0, 1, . . .

Proof.

Using the form of Ai (m) (x) given and using (62) we note that Ai (m+1) = p m Ai + q m Ai (x) + p m Ai + q m Ai = p m Ai + q m Ai (x) + p m Ai + xq m Ai, where we used the Airy differential equation Ai = xAi. This gives the recursion (73). The statement about the degree is easily checked for m = 0, 1. Suppose it holds upto a general m ≥ 2. Then

|p m | ∞ + |q m | ∞ ≤ Π m-1 =0 (d + 1) ≤ Π m+1 =2 d ≤ (d m+1 !) 2 .
Finally, we have for a polynomial p of degree d, and |x| < 1,

|p(x)| |p| ∞ ≤ 1 + |x| + • • • + |x| d ≤ 1 1 -|x| .
The last statement of the lemma is known for p = 0. Suppose it holds for p and use the Airy differential equation Ai = xAi. That gives Ai (3p+2) = xAi (3p) + (3p -1)Ai (3(p-1)+2) , which shows the claim.

Remark 9.1. The first few polynomials p m and q m in the theorem are given by Ai (2) (x) = xAi(x),

Ai (3) (x) = Ai(x) + xAi (x), Ai (4) (x) = x 2 Ai(x) + 2Ai (x), Ai (5) (x) = 4xAi(x) + x 2 Ai (x), Ai (6) (x) = (x 3 + 4)Ai(x) + 6xAi (x), Ai (7) (x) = 9x 2 Ai(x) + (x 3 + 10)Ai (x), Ai (8) (x) = (x 4 + 28x)Ai(x) + 12x 2 Ai (x).

Therefore, if p = 0, . . . , 3 and j 1 + • • • + j p = p ≤ 4p,

Z (p) (G ε (θ)) d j1 G ε (θ) dθ j1 • • • d jp G ε (θ) dθ jp ≤ C (1 + |εθ|) p 1 + ε p |θ| 4p (C ε) j1+•••+jp |θ| 4-j1+4-j2+•••+4-jp (1 + |εθ|) j1+•••+jp = C ε p |θ| 4p-p 1 + ε p |θ| 4p = C ε p /4 (ε 1/4 |θ|) 4p-p 1 + (ε 1/4 |θ|) 4p ≤ C ε p /4 .
From these estimates we get, with p = p = 0, . . . , 3,

|σ ε (θ)| = |Z(G ε (θ))| ≤ C, d dθ σ ε (θ) = Z (1) (G ε (θ)) dG ε (θ) dθ ≤ Cε 1/4 , d 2 dθ 2 σ ε (θ) = Z (2) (G ε (θ)) dG ε (θ) dθ 2 + Z (1) (G ε (θ)) d 2 G ε (θ) dθ 2 ≤ Cε 1/2 , d 3 dθ 3 σ ε (θ) = Z (3) (G ε (θ)) dG ε (θ) dθ 3 + 3Z (2) (G ε (θ)) d 2 G ε (θ) dθ 2
dG ε (θ) dθ + 3Z (1) (G ε (θ))

d 3 G ε (θ) dθ 3 ≤ Cε 3/4 .
This shows the lemma.

Figure 1 .

 1 Figure 1. The fold caustic: Example of solution with ray tracing picture for 0 ≤ x ≤ 1.

Figure 2 .

 2 Figure 2. A Gaussian beam starting at (x, y) = (0, -0.5) with direction η 0 = sin Θ = 3/4. The central ray γ is indicated with a solid line.

Figure 3 .

 3 Figure 3. The simplified map s * (y).

Thus, since φ

  a has no stationary points on K we can use Lemma B.1 with a(ζ) = b(ζ)ζ p and D an open bounded set containing K. This gives

Lemma 6 . 3 .

 63 Let I = k 1/3 ψ(θ/R)θ p e ikφa(δ,θ,η) dθ,where c 0 ≤ R with c 0 as in Lemma 5.3 and |δ| ≤ 1. If c 0 (1 + |η| 1/2 ) ≤ R and n ≥ 1 + p/2, there is a constant C n , independent of k, R, δ and η, such that

2 3 ikx 3/ 2 ⇒

 22 Ai(αk 2/3 X) = Ai(αk 2/3 x c )

  (62) Ai (z) = zAi(z), Bi (z) = zBi(z).(P2) Ai and Ai only have zeros on the negative real line. The zeros do not coincide. Ai(s) is positive and decreasing for s ≥ 0. (P3) Bi and Bi also only have zeros on the negative real line. The zeros do not coincide. Bi(s) is positive and increasing for s ≥ 0

Lemma 9 . 1 .

 91 Let β = e iθ with |θ| ≤ π/3. There is a constant C such that Ai(βs) -Ai(βs) ≤ Cs -3/2 Ai(βs) , s > 0. (66) Moreover, for each s 0 > 0 there is a constant C(s 0 ) such thatAi(βs) -Ai(βs) ≤ C(s 0 )s -3/2 |Ai(βs)| , s ≥ s 0 . (67) Moreover, for s ∈ R and α = exp(iπ/3) there is a constant C such that |Ai(s)| ≤ C(1 + |s|) -1/4 , (68) |Ai (s)| ≤ C(1 + |s|) 1/4 , (69) |Ai(αs)| ≥ C(1 + |s|) -1/4 , (70) |Ai (s) + i √ -sAi(s)| ≥ C (1 + |s|) 1/2 |Ai(s)|, s ≥ 0, (1 + |s|) 1/4 , s < 0. (71)Proof. We define Φ 0 by the relation Ai(z) = Ai(z)Φ 0 (z). It has an asymptotic expansion,Φ 0 (z) ∞ n=0 c n z -3n/2 , c n = (-1) n Γ(n + 5/6)Γ(n + 1/6)From the estimates on P and Q of [8, Appendix A,Lemma 7] one obtains the uniform estimate|Φ 0 (z) -1| ≤ |c 1 | |z| 3/2 ,valid for | arg(z)| ≤ π 3 . Then (66) follows directly with C = |c 1 |, Ai(βs) -Ai(βs) = Ai(βs) |Φ 0 (βs) -1| ≤ |c 1 |

2 , 2 ≤-1 4 .Ais 3 / 2

 22432 s 1 := (2|c 1 |) 2/3 we have |Φ 0 (sβ)| ≥ 1/2, and since Ai = AiΦ 0 only has zeros on the negative real line, there is a positive infinum, when s ≥ s 0 , proving (67) with C(s 0 ) = |c 1 |/d(s 0 ). Next, to show (68) and (69) we note first that, for real s > 0, Cs Then (64, 65) readily give |Ai(s)| ≤ Ai(s) + C 1 + Ai(s) |s| -3/2 ≤ C|s| -1/4 and Ai (s) ≤ |s| Ai(s) + C 1 + Ai(s) |s| -1 ≤ C|s| 1/4 . which extend to (68) and (69) as Ai(0) is bounded. The lower bound (70) follows for s ≥ s 0 from the previous estimates, Ai(αs) ≤ Ai(αs) -Ai(αs) + |Ai(αs)| ≤ |c 1 | Ai(αs) + |Ai(αs)| ,

1 .Lemma 9 . 2 .

 192 2.11] (72) Ai(αs) = ᾱ 2 [Ai(-s) + iBi(-s)], and (P3) implies that |Ai(αs)| ≥ |Bi(-s)|/2 = Bi(-s)/2 ≥ Bi(0)/2 > 0 for s ≤ 0. This gives the bound (70) also for s ≤ s 0 .For (71) we consider s ≥ 0 and use (64, 65)2 √ s| Ai(-s)| ≤ |2 √ s Ai(-s) -iAi (-s) -√ sAi(-s)| + |iAi (-s) + √ sAi(-s)| ≤ |2 √ s Ai(-s) -Ai (-s)| + √ s|2 Ai(-s) -Ai(-s)| + |iAi (-s) + √ sAi(-s)| ≤ Cs -1 + |iAi (-s) + √ sAi(-s)| = Cs -1 + |Ai (-s) + i √ sAi(-s)|. Hence, |Ai (-s) + i √ sAi(-s)| ≥ 2 √ s| Ai(-s)| -Cs -1 = Cs 1/4 -Cs -Since the zeros of Ai and Ai do not coincide (no double roots) we get|Ai (-s) + i √ sAi(-s)| = |Ai (-s)| + √ s|Ai(-s)| = 0, s > 0,and the estimate (71) for s < 0 follows. Moreover, by (P2), when s ≥ 0,|Ai (s) + i √ -sAi(s)| = |Ai (s) -√ sAi(s)| = √ sAi(s) -Ai (s) = √ s|Ai(s)| + |Ai (s)| ≥ max √ s|Ai(s)|, |Ai(s)| min 0≤t≤s |Ai (t)| |Ai(t)| ≥ 1 2 ( √ s + C)|Ai(s)|,which gives (71) for s ≥ 0. Here we also used the fact that lim s→+∞|Ai (s)| |Ai(s)| = lim s→+∞ | √ s Ai(s)| | Ai(s)| = ∞ by (64) and (65). For the Airy function we have Ai (m) (x) = p m (x)Ai(x) + q m (x)Ai (x),

  The degree of their sum satisfies deg(p m + q m ) = m/2 and, for |x| < 1,|p m (x)| + |q m (x)| ≤ (d m+1 !)

  deg(p m+1 + q m+1 ) = deg(p m + q m + p m + xq m ) = deg(p m + xq m ) = deg(p m-1 + xq m-1 + xp m-1 + xq m-1 ) = deg(x(p m-1 + q m-1 )) = deg(p m-1 + q m-1 ) + 1.That deg(p m + q m ) = d m follows by induction. Since the polynomials all have positive coefficients, it also follows thatd m = max(deg(p m ), deg(q m )). For a polynomial p, let |p| ∞ denote its largest coefficient in magnitude. Then |xp| ∞ = |p| ∞ and |p | ∞ ≤ deg(p)|p| ∞ . Consequently, |p m+1 | ∞ + |q m+1 | ∞ = |p m + xq m | ∞ + |p m + q m | ∞ ≤ |p m | ∞ + |xq m | ∞ + |p m | ∞ + |q m | ∞ ≤ deg(p m )|p m | ∞ + |q m | ∞ + |p m | ∞ + deg(q m )|q m | ∞ ≤ (max(deg(p m ), deg(q m )) + 1)(|p m | ∞ + |q m | ∞ ) = (d m + 1)(|p m | ∞ + |q m | ∞ ). Therefore, since |p 0 | ∞ + |q 0 | ∞ = 1,

3

  ) 1/4 ≤ C 3 6 1/4 k 5/6 η 2 . It follows that |R 1 | ≤ M kη 2 where M = 2 √ π6 1/4 C 3 with C 3 beingthe constant in (69).There is a constant M dependent on x, but independent of η and k ≥ 1, such that

	8.2. Estimate of R 2 .
	Lemma 8.4.

  LAGA, UMR7539 Université Sorbonne Paris Nord, 99, Avenue J. B. Clement, 93430 Villetaneuse, France IRL CNRS-CRM, Université de Montréal, Chemin de la tour, Montréal, Canada. Department of Mathematics, KTH, 10044 Stockholm, Sweden

By appealing to Lemma C.1 with ε = k -1/3 we conclude that ||z(•/k 1 3 , k)|| W 3,∞ (R) is bounded uniformly for k ≥ 1. Consequently, we can use Lemma 6.2 with n = (12 + 1)/5 = 3 together with (37) and (38) to show that

We have thus proved that (58)

|I main -I S | ≤ Ck -1 .

From Lemma 6.3 we know that I p ≈ 2πAi (p) (ρ)i p /k p/3 and we therefore introduce the approximation ĨS of I S obtained by replacing I p with the corresponding Airy function,

. By (56) we can use Lemma 6.3 with large enough n to obtain

The next step is to show that ĨS is close to Ai(ρ). Upon using the identities for Ai (m) with m = 2, 4, 5, 8 given in Remark 9.1 we can simplify the expression for ĨS as follows

From ( 55), ( 7) and ( 16), we obtain

Since ρ is bounded by (53) and Ai is smooth around ρ = 0, this shows that

Note that the dependence on k 2 η 2 which appears here, also appears in the estimate of R 2 in Lemma 8.4.

Appendix A. Proof of Lemma 2.1

We show this for k = 1 so that F k = F. The case with general k follows from a simple rescaling. Let φ ∈ S be a test function and • , • the duality pairing between S and S . Then since f ψ t ∈ L ∞ ⊂ S and, by dominated convergence,

where we used the facts that F(φ) ∈ S ⊂ L 1 , |ψ t | ≤ 1 for all t and f ψ t → f pointwise. This is true for all φ ∈ S and therefore F(f ψ t ) → F(f ) in S , proving the first statement.

That F((f * g)ψ t ) → F(f )F(g) follows from the first statement since f * g ∈ L ∞ when g ∈ S and F(f * g) = F(f )F(g). The last part of the second statement is true, since F(g)φ ∈ S and therefore the first part gives

This shows the lemma.

Appendix B. Proof of the non-stationary phase identities

Below is a proof of identities used in the non-stationary phase lemma. The identities show how the rewritten integral depends on the derivatives of the phase function. In order to do that we use the spaces of functions defined in (35) and (36). We thus need to show that there exist u ,n such that

When n = 0 this simply says that u 0,0 = 1 ∈ U 0 (φ). Suppose the claim holds for n and consider

For the first term we have

For the third term

For the second term, consider one basis function w ∈ W p (φ),

for some M . Then

Hence, u ,n ∈ U n+1-(φ), as it is a linear combination of derivatives of functions in U n-(φ). This shows that L n+1 [a] is of the correct form and the lemma is proved.

Appendix C. Boundedness of Z

Here we consider the scaled remainder term in the Taylor expansion of exp(iz), ( 76)

We have the following lemma.

Then there is a constant C such that

Proof. We begin by estimating Z and its first three derivatives for z with non-negative imaginary part. Then all derivatives of Z are bounded, since

Furthermore, from the first part of the definition (76) we get, for p = 0, . . . , 3 and |z| ≥ 1, (3,p+1) .

It follows that there is a constant C such that