General superpositions of Gaussian beams and propagation errors

TL;DR

This paper extends Gaussian beam superpositions to arbitrary phase space sets, deriving propagation error estimates and analyzing their sharpness, thereby broadening the applicability of high-frequency PDE solutions.

Contribution

It introduces a superposition of Gaussian beams over arbitrary phase space sets and applies recent tools to estimate propagation errors, advancing high-frequency solution methods.

Findings

Propagation error of order $k^{1- rac{N}{2}- rac{d-m}{4}}$ derived

Error estimate sharpness analyzed through examples

Method applicable to general high frequency PDE solutions

Abstract

Gaussian beams are asymptotically valid high frequency solutions concentrated on a single curve through the physical domain, and superposition of Gaussian beams provides a powerful tool to generate more general high frequency solutions to PDEs. We present a superposition of Gaussian beams over an arbitrary bounded set of dimension $m$ in phase space, and show that the tools recently developed in [ H. Liu, O. Runborg, and N. M. Tanushev, Math. Comp., 82: 919--952, 2013] can be applied to obtain the propagation error of order $k^{1- \frac{N}{2}- \frac{d-m}{4}}$, where $N$ is the order of beams and $d$ is the spatial dimension. Moreover, we study the sharpness of this estimate in examples.

Tables

Table 1. Table 1. 3D Gaussian beam single point errors and orders of convergence.

$t$	k=320	k=640		k=1280		k=2560
	error	error	order	error	order	error	order
0.4	0.109724	0.064004	0.78	0.0347248	0.88	0.0177462	0.97
0.55	0.0820207	0.0420894	0.96	0.0213659	0.98	0.0118253	0.85
0.7	0.0822853	0.0418797	0.97	0.0211195	0.99	0.0102407	1.04

Table 2. Table 2. 3D Gaussian beam energy errors and orders of convergence.

$t$	k=320	k=640		k=1280		k=2560
	error	error	order	error	order	error	order
0.4	0.111507	0.0671302	0.73	0.0354353	0.92	0.0185048	0.94
0.5	0.0716308	0.0388652	0.88	0.0193636	1.01	0.00994688	0.96
0.55	0.0825064	0.0429692	0.94	0.0213441	1.01	0.0108103	0.98
0.7	0.0834459	0.0427814	0.96	0.0211234	1.02	0.0106206	0.99
0.8	0.0458945	0.0242241	0.92	0.0100285	1.27	0.0053213	0.91

Table 3. Table 3. L ∞ superscript 𝐿 L^{\infty} errors and orders of convergence of 2D Gaussian beam superposition.

$t$	k=80	k=160		k=320
	error	error	order	error	order
0.15	0.0134258	0.00775526	0.79	0.00390146	0.99
0.30	0.0349802	0.0205824	0.77	0.0119712	0.78
0.40	0.0439910	0.0253134	0.80	0.0155627	0.70
0.42	0.0469892	0.0241403	0.96	0.0155593	0.63
0.45	0.0514442	0.0258481	0.99	0.0147372	0.81
0.50	0.0578684	0.0296832	0.96	0.0161166	0.88
0.60	0.0694715	0.0371298	0.90	0.0184097	1.01
0.70	0.0795142	0.0437080	0.86	0.0229725	0.93
0.80	0.0878869	0.0485395	0.86	0.0254761	0.93

Table 4. Table 4. Energy errors and orders of convergence of 2D Gaussian beam superposition.

$t$	k=80	k=160		k=320
	error	error	order	error	order
0.15	0.0138	0.0071	0.96	0.0034	1.06
0.30	0.0274	0.0141	0.96	0.0068	1.05
0.40	0.0364	0.0187	0.96	0.009	1.06
0.42	0.0382	0.0196	0.96	0.0094	1.06
0.45	0.0408	0.021	0.96	0.0101	1.06
0.50	0.0452	0.0233	0.96	0.0112	1.06
0.60	0.054	0.0278	0.96	0.0134	1.05
0.70	0.0626	0.0323	0.95	0.0156	1.05
0.80	0.071	0.0367	0.95	0.0177	1.05

Equations

$$Pu={\partial }_t^{2}u(x,t)-c(x{)}^2\Delta u(x,t)=0,(x,t)\in {\mathbb{R}}_x^d\times {\mathbb{R}}_t$$

$$(u(x,0),{\partial }_{t}u(x,0))=(B_0(x),k{B}_1(x))e^{ik{S}_0(x)},$$

$$v(x,t)=(a_0(x,t)+k^{-1}{a}_1(x,t)+\cdots +k^{-p}{a}_p(x,t))e^{ik\varphi (x,t)}$$

$$[Pv](x,t)=e^{ik\varphi (x,t)}\sum _{j=0}^{p+1}{k}^{2-j}{L}_j(x,t).$$

$$u_{GB}(x,t)=k^{\frac{m}{2}}{\int }_{K_0}v(x,t;X_0)d{X}_0,$$

$$||u|{|}_E^2=\frac{1}{2}{\int }_{{\mathbb{R}}^n}{c}^{-2}(x)|{\partial }_{t}u{|}^2+|{\nabla }_{x}u{|}^{2}dx$$

$$||u(\cdot ,t)-u_{GB}(\cdot ,t)|{|}_E\le C(T)k^{1-N/2}$$

$$||u(\cdot ,t)-u_{GB}(\cdot ,t)|{|}_E\le C(T)k^{1-N/2-(d-m)/4}$$

$$||u(\cdot ,t)-u_{GB}(\cdot ,t)|{|}_E\le C(T)k^0$$

$$Pu=0,$$

$$\Vert (u-v)(\cdot ,t_2){\Vert }_S\le \Vert (u-v)(\cdot ,t_1){\Vert }_S+C{k}^q{\int }_{t_1}^{t_2}\Vert Pv(\cdot ,\tau ){\Vert }_{L^2}d\tau .$$

$$\Vert u(\cdot ,t){\Vert }_E={\left(\frac{1}{2}{\int }_{{\mathbb{R}}^d}\left(c(x{)}^{-2}|{\partial }_{t}u(x,t){|}^2+|{\nabla }_{x}u(x,t){|}^2\right)dx\right)}^{1/2},$$

$$P=-iϵ{\partial }_t+\frac{{ϵ}^2}{2}\Delta$$

$$\Vert u(\cdot ,t_2){\Vert }_S=\Vert u(\cdot ,t_1){\Vert }_S,\forall t_1,t_2.$$

$$C{k}^{-\alpha }\ge \Vert (u-u_{GB})(\cdot ,t_1){\Vert }_S\text{and}\Vert (u-u_{GB})(\cdot ,t_2){\Vert }_S\ge c{k}^{-\beta },$$

$$\Vert (w_1-u_{GB})(\cdot ,t_1){\Vert }_S=0\; and\; \Vert (w_1-u_{GB})(\cdot ,t_2){\Vert }_S\ge c_0{k}^{-\beta }$$

$$\Vert (w_2-u_{GB})(\cdot ,t_1){\Vert }_S\ge c_0{k}^{-\beta }and\; \Vert (w_2-u_{GB})(\cdot ,t_2){\Vert }_S=0$$

$$\Vert (u-w_1)(\cdot ,t_2){\Vert }_S=\Vert (u-w_1)(x,t_1){\Vert }_S=\Vert (u-u_{GB})(\cdot ,t_1){\Vert }_S\le C{k}^{-\alpha }.$$

$$\Vert u_{GB}(\cdot ,t_2)-w_1(\cdot ,t_2){\Vert }_S$$

$$\ge c{k}^{-\beta }-C{k}^{-\alpha }.$$

$$\Vert (u-w_2)(\cdot ,t_1){\Vert }_S=\Vert (u-w_2)(x,t_2){\Vert }_S=\Vert (u-u_{GB})(\cdot ,t_2){\Vert }_S\ge c{k}^{-\beta }.$$

$$\Vert u_{GB}(\cdot ,t_1)-w_2(\cdot ,t_1){\Vert }_S$$

$$\ge c{k}^{-\beta }-C{k}^{-\alpha }.$$

$$u_{GB}(x,t)=k^{m/2}{\int }_{K_0}v(x,t;X_0)d{X}_0,$$

$$v(x,t;X_0)=\sum _{j=0}^{\lceil N/2\rceil -1}{k}^{-j}{\rho }_{\eta }(x-x(t;X_0))a_j(t,x-x(t;X_0))e^{ik\varphi (t,x-x(t;X_0))},$$

$$\rho_{\eta}(z)=\left\{\begin{array}[]{ll}1&|z|\leq\eta,\\ 0&|z|\geq 2\eta\end{array}\right.\quad 0<\eta<\infty.$$

$$\varphi (t,y)=S(t;X_0)+p(t;X_0)\cdot y+\frac{1}{2}y\cdot M(t;X_0)y,$$

$$\dot{x}$$

$$\dot{S}$$

$$\dot{M}$$

$$\Vert u_{GB}-u{\Vert }_S\lesssim k^{1-N/2-(d-m)/4},$$

$$\Vert u_{GB}(\cdot ,t)-u(\cdot ,t){\Vert }_S\lesssim \Vert u_{GB}(\cdot ,0)-u(\cdot ,0){\Vert }_S+k^{1-N/2-(d-m)/4}.$$

Full text

1. Introduction

In this paper we investigate issues related to the accuracy of Gaussian beam approximations to high frequency wave propagation. This is related to recent results on Gaussian beam methods in [5, 6, 7, 8, 9, 10, 11, 13]. Our model equation is the acoustic wave equation

[TABLE]

where $c(x)$ is a positive smooth function. The initial data are given by

[TABLE]

where $k\gg 1$ and $\nabla S_{0}\neq 0$, so that the data are highly oscillatory. Propagation of high frequency oscillations leads to mathematical and numerical challenges in solving wave propagation problems.

We study the errors which arise when one approximates solutions to the initial value problem (1.1) by superpositions of Gaussian beams. Our starting point is [9], and we refer the reader to it for more references to earlier results on superpositions of beams. In addition, some recent effort has also been made to extend the Gaussian beam method to more complex settings such as symmetric hyperbolic systems with polarized waves [3], the Schrödinger equation with discontinuous potentials [4], and wave equations in bounded convex domains [1, 2].

To compare the results in [9] with what we do here we need to recall some conventions. For a Gaussian beam

[TABLE]

we say that $v$ is an $N$th-order approximation to a solution of $Pu=0$ when the sequence of equations (from geometric optics) $L_{j}(x,t)=0$, $j=0,1,\dots$, holds to order $N+2-2j$ on the central ray path, where

[TABLE]

Analogously to [9] we use superpositions of the form

[TABLE]

where $K_{0}$ is a submanifold of dimension $m$ in phase space that does not intersect $\{(x,p):\;p=0\}$, and the central ray for $v(x,t;X_{0})$ has initial data $(x(0),p(0))=X_{0}$. In this paper we are considering superpositions over submanifolds of $2d$-phase space of dimension at most $d$. Finally we use the unscaled energy norm

[TABLE]

in place of the scaled energy norm in [9] which has an additional factor of $k^{-1}$. With these conventions the principal result of [9] becomes

Theorem 1.1**.**

([9]) When $u(x,t)$ is the exact solution to $Pu=0$ with the initial data of the superposition $u_{GB}$ of Gaussian beams of order $N$ over a compact subset $K_{0}$ of dimension $d$ in $\mathbb{R}^{d}$, the error estimate

[TABLE]

holds for $t\in[0,T]$.

In this note we extend that to

Theorem 1.2**.**

With the hypotheses in Theorem 1.1

[TABLE]

when $K_{0}$ is a bounded domain in phase space of dimension $m$.

Comparing Theorems 1.1 and 1.2 one sees that Theorem 1.1 is the special case where $K_{0}$ is a domain in $\mathbb{R}^{d}$ and hence $m=d$. In Theorem 1.2 the initial data is not restricted to the “WKB” form in (1.2). In this paper we will always use superpositions of the form (1.3) with beams that have leading amplitudes independent of $k$. Later in this paper we sometimes fix the dependence of the error on $k$ by dividing by the energy norm of the initial data. The decrease in the error becomes faster as $m$ decreases. This might be counter-intuitive, but it is consistent with the results in §5 of [9] where for a single first order beam ($N=1$ and $m=0$) in 2 dimensions

[TABLE]

in the unscaled energy norm above.

Theorem 1.2 is sharp in some cases. In Section 4 we give an example with $d=3$, $m=2$ and $N=1$, where the error as a function of $k$ decays no faster than the rate in (1.5). However, the initial data in this example is not of the form (1.2). A question that was left open in [9] is whether (1.4) is sharp for data of that form. Numerical evidence in [9] suggests that it is sharp when $N$ is even, but that when $N$ is odd the exponent on $k$ should be decreased by $1/2$, giving a faster decrease in the error as $k$ increases. There are partial results on this conjecture. For superpositions of first order beams (N=1) for the semi-classical Schrödinger equation a proof of the faster decay of the error in $L^{2}$ is presented in [13]***Zheng’s method can be applied to estimate errors in gaussian beam approximations for the acoustic wave equation in the $L^{2}$ and energy norms. This is consistent with the results in Table 2 and 4., based on ideas from [12].

For both the wave equation and the semi-classical Schrödinger equation, in [11] the authors show that, away from caustics, the error has, uniformly, the faster decay rate in the maximum norm. However, close to caustics, their estimate degenerates.

This paper is organized as follows: In Section 2 we derive a lower bound on the error for approximation by beam superpositions using energy conservation. In Section 3 we prove Theorem 1.2. In Section 4 we construct the example mentioned above. In Section 5 we construct a superposition with $N=1$ for the acoustic wave equation with initial data of the form (1.2) that develops a focus caustic at the origin. In a numerical study of this example we see the faster decay in the error conjectured in [9]. In Section 6 we construct an example in two space dimensions which develops a fold caustic on the unit circle. Here we again see numerically the faster decay conjectured in [9] in the energy norm, but in the maximum norm the decay is slower at some times. Section 7 is concerned with initial asymptotic rates shown by the construction of various examples. These examples illustrate the initial data that can arise from superpositions of the form (1.3) and their respective energy norms. Some final remarks are given in section 8. Notation: Throughout this paper, we use the notation $A\lesssim B$ to indicate that $A$ can be bounded by $B$ multiplied by a constant independent of the frequency parameter $k$. $A\sim B$ stands for $A\lesssim B$ and $B\lesssim A$.

2. Energy conservation and lower error bound

Error estimates for Gaussian beam superpositions are based on the well-posedness of the underlying equation. For an equation of the form

[TABLE]

we recall the well-known results here (see, e.g., [9])

Theorem 2.1**.**

Let $u$ be an exact solution of the wave equation (2.1), and $v$ be an approximate solution of the same problem, then we have the generic well-posedness estimate

[TABLE]

These apply to both

• •

the wave equation with $P=\partial_{t}^{2}-c^{2}(x)\Delta$, $q=0$, and $\|\cdot\|_{S}$ is the energy norm

[TABLE]

• •

and the Schrödinger equation with $q=1$, $\epsilon=\frac{1}{k}$,

[TABLE]

and $\|\cdot\|_{S}$ is the standard $L^{2}$ norm.

The lower bound on approximation errors is a consequence of the conservation law

[TABLE]

Theorem 2.2**.**

Let $u_{GB}$ be a Gaussian beam superposition, and $u$ be an exact solution of $Pu=0$. Assume that for some $\alpha>\beta>0$ there are times $t_{1}$ and $t_{2}$ and positive constants $C,\ c$ such that for $k\geq 1$

[TABLE]

then there are exact solutions $w_{1}$ and $w_{2}$ and a $c_{0}>0$ such that

[TABLE]

and

[TABLE]

for $k$ sufficiently large.

Proof.

Let $w_{1}(x,t)$ be the exact solution with data at $t=t_{1}$ that agree with the data of $u_{GB}(x,t)$ at $t=t_{1}$. By (2.3) we have

[TABLE]

It follows that

[TABLE]

For the other case, we argue the following manner. Let $w_{2}(x,t)$ be the exact solution with data at $t=t_{2}$ that agree with the data of $u_{GB}(x,t)$ at $t=t_{2}$. By energy conservation we have

[TABLE]

It follows that

[TABLE]

∎

Remark 2.1*.*

This result may be used to identify the source of accuracy loss of the Gaussian beam superposition or other types of approximate solutions.

3. Propagation error of Gaussian beam superpositions

Let $K_{0}$ be an arbitrary bounded set in phase space with dimension $m$. Given a point $X_{0}\in K_{0}$, we denote the $N$-th order Gaussian beam as $v(x,t;X_{0})$, if we let $X_{0}$ range over $K_{0}$, we can form a superposition of Gaussian beams,

[TABLE]

as an approximation to the exact solution for wave equation (2.1) with initial data $u_{GB}(x,0)$.

We recall that the general form of the N-th order Gaussian beam defined in [9] is

[TABLE]

where $\rho_{\eta}(\cdot)\geq 0$ is a smooth cutoff function satisfying $\rho_{\infty}=1$ and

[TABLE]

In this construction the parameter $\eta$ is chosen as $\eta=\infty$ for the first order superposition and it is taken small enough to make $Im(\phi(t,y))\geq\delta|y|^{2}$ for $t\in[0,T]$ and $|y|\leq 2\eta$ for higher order superpositions.

For first order beams,

[TABLE]

associated with the first several ODEs defined by

[TABLE]

For equation (1.1), $H(x,p)=\pm c(x)|p|$, for which two wave modes need to be included in the superposition. We assume that $K_{0}$ does not intersect $\{(x,p)|\;p=0\}$. No such assumption is needed for the Schrödinger equation with $H(x,p)=\frac{1}{2}|p|^{2}$. These construction details will not be used in our error analysis, but may be helpful as a reference for reading examples constructed in sections 4-6.

We now state the main result of the propagation error for superposition (3.1).

Theorem 3.1**.**

Let $u_{GB}$ be the Gaussian beam superposition defined in (3.1) based on $N$-th order beams emanating from a compact subset of the $m$-dimensional manifold $K_{0}$ in phase space, and $u$ be the exact solution to $Pu=0$ subject to the initial data $u_{GB}(x,0)$, we then have the following estimate on the propagation error,

[TABLE]

where $m$ is the dimension of the domain on which initial beams are sampled, and $d$ is the spatial dimension.

Remark 3.1*.*

Note that operator $P$ is initially defined in (1.1), but also used for Schrödinger operator in Section 2. This theorem includes the proof of Theorem 1.2, but it is also valid for the Schrödinger equation due to the basic estimate (2.2) and the estimate of $\|Pu_{\rm GB}\|$ to be carried out in this section.

We proceed to complete the proof of this theorem by following the general steps as in the proof of [9, Theorem 1.1]. The main difference here is that the initial set $K_{0}$ can be rather arbitrary in phase space. The way that distance between beams is measured must here be allowed to vary smoothly with the beam’s initial point in phase space.

Before we outline the proof of the above result, we present a result, which shows that the accuracy of the initial approximation can be treated separately.

Corollary 3.2**.**

Let $u_{GB}$ be the Gaussian beam superposition defined in (3.1) based on $N$-th order beams, and $u$ be the exact solution to $Pu=0$ subject to a given initial data $u(x,0)$, then

[TABLE]

Proof.

Let $w$ be another exact solution with initial data $u_{GB}(x,0)$, then we have

[TABLE]

The energy conservation tells that

[TABLE]

These combined with the triangle inequality

[TABLE]

lead to (3.3). ∎

In this section, we focus only on the residual error, where the residual can be written (following the notation of Liu, Runborg and Tenushev [9] and Liu, Ralston, Runborg and Tanushev[10]) in the form

[TABLE]

where $Pv(x,t;X_{0})$ is a finite sum of terms of the form

[TABLE]

with bounds

[TABLE]

Here $g$ is smooth and supported or at least bounded on

[TABLE]

and $\phi$ is the $N$-th order Gaussian beam phase. Here $\tilde{\eta}$ is chosen as a small number for first order beams, but can be taken as $\eta$ for higher order beams. Moreover, $O(k^{-\infty})$ indicates terms exponentially small in $1/k$. After neglecting these terms and using (3.4) we can bound the $L^{2}$ norm of $P[u_{GB}]$ by

[TABLE]

where the term $I$ is of the form

[TABLE]

Here

[TABLE]

The function $g$ and its derivatives are bounded, for $0\leq t\leq T$,

[TABLE]

The rest of this section is dedicated to establishing the following inequality

[TABLE]

With this estimate we have

[TABLE]

which together with the well-posedness estimate (2.2) leads to the desired estimate (3.2).

Lemma 3.3** (Non-squeezing lemma).**

Let $X=(x(t;X_{0}),p(t;X_{0}))$ be the Hamiltonian trajectory starting from $X_{0}\in K_{0}$ with $K_{0}$ bounded. Assume that $X(0;X_{0})\in C^{2}(K_{0})$. Then,

[TABLE]

The non-squeezing lemma [9] says that the distance in phase space between two smooth Hamiltonian trajectories will not shrink from its initial distance. Here one may take any $l^{p}$ distance since from $X-X^{\prime}=(x-x^{\prime},0)+(0,p-p^{\prime})$ we have

[TABLE]

We recall some main estimates from [9] for proving (3.7).

Lemma 3.4** (Phase estimates).**

Let $\tilde{\eta}$ be small and $x\in D(\tilde{\eta},X_{0},X_{0}^{\prime})$ with

[TABLE]

• •

For all $X_{0},X_{0}^{\prime}\in K_{0}$ and sufficiently small $\tilde{\eta}$, there exists a constant $\delta$ independent of $k$ such that

[TABLE]

• •

For $|\gamma(x,t;X_{0})-\gamma(x,t;X_{0}^{\prime})|<\theta|X_{0}-X_{0}^{\prime}|$,

[TABLE]

where $C(\theta,\tilde{\eta})$ is independent of $x$ and positive if $\theta$ and $\tilde{\eta}$ are sufficiently small.

Decompose $I$ as

[TABLE]

with

[TABLE]

where $\chi_{j}(x,t,X_{0},X_{0}^{\prime})\in C^{\infty}$ is a partition of unity such that

[TABLE]

We first estimate $I_{1}$, which corresponds to the non-caustic region of the solution.

[TABLE]

Here we have used the fact that $|\gamma-\gamma^{\prime}|>\theta|X_{0}-X_{0}^{\prime}|$ on the support of $\chi_{1}$. For the inner integral over $D=\Omega(\tilde{\eta};X_{0})\cap\Omega(\tilde{\eta};X_{0}^{\prime})$, we have

[TABLE]

From this it follows that

[TABLE]

Letting $\Lambda=\sup_{X_{0},X_{0}^{\prime}\in K_{0}}|X_{0}-X_{0}^{\prime}|<\infty$ be the diameter of $X_{0}$, we have

[TABLE]

which concludes the estimate of ${\mathcal{I}}_{1}$.

In order to estimate ${\mathcal{I}}_{2}$ we use a version of the non-stationary phase lemma.

Lemma 3.5** (Non-stationary phase lemma).**

Suppose that $u(x;\zeta)\in C_{0}^{\infty}(\Omega\times Z)$, where $\Omega$ and $Z$ are compact sets and ${\psi(x;\zeta)\in C^{\infty}(O)}$ for some open neighborhood $O$ of $\Omega\times Z$. If $\nabla_{x}\psi$ never vanishes in $O$, then for any $K=0,1,\ldots$,

[TABLE]

where $C_{K}$ is a constant independent of $\zeta$.

We now define

[TABLE]

Non-stationary phase Lemma 3.5 can be applied to ${\mathcal{I}}_{2}$ with $\zeta=(X_{0},X_{0}^{\prime})\in K_{0}\times K_{0}$ to give,

[TABLE]

On the support of $\chi_{2}$ the difference $|X_{0}-X_{0}^{\prime}|$ can be arbitrary small, in which case this estimate is not useful. Following [9], we use the fact that the estimate is true also for $K=0$ so that ${\mathcal{I}}_{2}$ can be bounded by the minimum of the $K=0$ and $K>0$ estimates. Therefore,

[TABLE]

Finally, letting $\Lambda=\sup_{X_{0},X_{0}^{\prime}\in K_{0}}|X_{0}-X_{0}^{\prime}|<\infty$ be the diameter of $K_{0}$, we compute

[TABLE]

if we take $K>m$. This shows the ${\mathcal{I}}_{2}$ estimate, which proves claim (3.7).

4. Example of a Gaussian beam superposition

Let $r=|x|$, $x\in\mathbb{R}^{3}$. Then for any smooth function $f$,

[TABLE]

satisfies $\partial_{t}^{2}u=\Delta u$. Take $f(r)=\exp(-ikr-kr^{2}/2)/k$. Then

[TABLE]

The exact solution here is a highly oscillatory spherical wave which concentrates on $r=|t|$ as $k\to\infty$. The Cauchy data of this solution at $t=0$ can be approximated very well by a superposition of Gaussian beams.

Note that

[TABLE]

since the integral is a radial solution of $\Delta w+k^{2}w=0$, which equals $4\pi$ at $x=0$, then we have

[TABLE]

where

[TABLE]

Let us approximate $u(x,t)$ by a superposition of beams

[TABLE]

Hence $u_{GB}(x,0)=u(x,0)$. It will turn out, somewhat surprisingly, that $\partial_{t}u_{GB}(x,0)$ is very close to $\partial_{t}u(x,0)$. In fact, the first order Gaussian beam can be explicitly given as

[TABLE]

where

[TABLE]

and $2\pi a(t)=i(1+it)^{-1}$. Note that $\partial_{t}v=(ik\partial_{t}\phi a+\partial_{t}a)e^{ik\phi},\hbox{ so }$

[TABLE]

Now we can compute

[TABLE]

Using

[TABLE]

we have

[TABLE]

Note that the first two terms in that expression equal $\partial_{t}u(x,0)$. To estimate the data we use the standard energy norm $||(u,\partial_{t}u)||_{E}^{2}=\int_{\mathbb{R}^{3}}|\partial_{t}u|^{2}+|\nabla_{x}u|^{2}dx$. In that norm the difference of the initial data satisfies

[TABLE]

So the relative error in the initial data is $O(k^{-3/2})$.

Now we get to the main point: How large is $u(x,t)-u_{GB}(x,t)$? We need to compute

[TABLE]

Introducing spherical coordinates so that $x\cdot\omega=|x|\cos\rho$ and $d\omega=\sin\rho d\rho d\phi$ with the domain of integration $0\leq\rho\leq\pi$ and $0\leq\phi\leq 2\pi$ and setting $|x|=r$, this becomes - after substituting $s=\cos\rho$

[TABLE]

where

[TABLE]

Note that, for $t>0$, the real part of the exponent in the integrand is strictly negative unless $s=1$ and $r=t$. Moreover, for $t>0$ and $r$ in a sufficiently small neighborhood of $t$ the maximum of the real part of exponent for $-1\leq s\leq 1$ is assumed at $s=1$. So we can find $u_{\rm GB}(x,t)$, up to terms of order $k^{-1}e^{-k(r-t)^{2}/2}$, by using the leading term in the integration by parts expansion: Choosing $\rho$ with support near $s=1$ and $\rho(1)=1$,

[TABLE]

One continues this expansion by repeated integration by parts. In particular, the integral term on the right is $O(k^{-2}e^{-k(r-t)^{2}/2})$. Since $\phi(\pm 1)=\pm r-t+i(r\mp t)^{2}/2$ and

[TABLE]

hence for $t>\delta>0$ and $r$ close to $t$,

[TABLE]

At this point we want to obtain a lower bound on $||u_{GB}(\cdot,t)-u(\cdot,t)||_{E}$. The dominant terms in the first derivatives of $u_{GB}(\cdot,t)-u(\cdot,t)$ come from the factor $\exp(ik(r-t))$ and bring down a factor of $k$. So, letting $s=r-t$, this leaves a dominant term which is a nonvanishing multiple of $se^{-ks^{2}/2}$, and hence has $L^{2}$ norm bounded below by a multiple of $k^{-3/4}$. That implies $||u_{GB}(\cdot,t)-u(\cdot,t)||_{E}\sim k^{-3/4}$. However, here the Gaussian beam superposition is missing a factor of $k$ compared to Theorem 1.2. Hence this example shows that Theorem 1.2 is sharp when $d=3$, $m=2$ and $N=1$.

5. An example for the 3D acoustic wave equation

This will be the construction of a Gaussian beam superposition for the initial value problem

[TABLE]

where $a(r)=0$ in a neighborhood of $r=0$. From here on $|x|=r$ will be used.

The exact solution to this initial value problem is

[TABLE]

extended to $\mathbb{R}$ by $f(s)=-f(-s)$. Note that for $t>0$

[TABLE]

and the solution has a strong peak at $r=0$, when $t$ is in the support of $a$. We want to see the effect of this caustic.

Following the “standard procedure” for first order beams, the Gaussian beam superposition will be

[TABLE]

where $\Phi^{+}(x,0;y)=\Phi^{-}(x,0;y)$, $\partial_{t}\Phi^{+}(x,0;y)=-\partial_{t}\Phi^{-}(x,0;y)$ and $A^{+}(0;y)=A^{-}(0;y)$. So we have two families of Gaussian beams

[TABLE]

where both phases $\Phi^{\pm}$ are based on the initial phase $S(x)=|x|$, but the $v^{\pm}$ are concentrated on the rays $(x(t),t)=(y\pm ty/|y|,t)$; see e.g., [7, superposition (3.1)].

From here on we will often use $y=s\omega$, $|\omega|=1$. Again the standard construction gives

[TABLE]

where $M(0;y)=(1/|y|)P_{\omega^{\perp}}+iI$ and $\partial_{t}M+MP_{\omega^{\perp}}M=0$. Here $I-P_{\omega^{\perp}}$ is the orthogonal projection on the span of $\omega$. A modest amount of computation shows

[TABLE]

where

[TABLE]

So

[TABLE]

The amplitudes $A^{\pm}$ are given by

[TABLE]

Since $x$ appears in $v^{\pm}$ only as $|x|$ and $x\cdot\omega$, we have $u_{GB}(x,t)=w(r,t)$. This can be seen by integrating in spherical coordinates. Also $\partial_{t}u_{GB}(x,0)=0$.

Now we need to determine the order of $||u(\cdot,t)-u_{GB}(\cdot,t)||_{E}$. The contributions to $u_{GB}$ from $\int A^{+}(t;y)\exp(ik\Phi^{+}(x,t;y))dy$ will be concentrated at $x=(t+s)\omega$, and, since $s\geq 0$ and we consider $t>0$, they will be negligible near $x=0$. Hence we will omit that term from all formulas from here on. Let $v={x\cdot\omega\over|x|}=\cos(\theta)$. While this is undefined at $x=0$, substitution of $v$ for $\theta$ in (5.2) leads to an integral in spherical coordinates that is well-behaved as $x\to 0$. Namely

[TABLE]

where $C=1-i(s-t)$ and $D=b(-t,s)/2$. Presumably one could evaluate this formula further, but that is a daunting calculation. Instead we offer the numerical results in the next section.

5.1. Numerical results

In this section and in the numerical results in Section 6.1 we will use relative norms to estimate errors, i.e. norms scaled by the corresponding norm of the beam superposition. In these examples that has the effect of decreasing the power of $k$ in the energy norm by one, and leaving the power unchanged in the $L^{2}$-norm, but Example 5 in section 7 shows that this does not always happen. Since the energy norm of the initial data is of order $k$ in both cases and $m=d$, Theorem 1.2 predicts a relative error of order $k^{-1/2}$ for first order beams. We will see that the actual error is numerically of order $k^{-1}$ as conjectured in [9].

We take $a(s)=4(s-r_{0})^{4}(s-r_{1})^{4}$ for $r_{0}\leq s\leq r_{1}$; $a(s)=0$ otherwise, here $r_{0}=0.1,\ r_{1}=1.0$. The evaluation of (5.2) is done using $80\times 80$ meshes of $[r_{0},r_{1}]\times[-1,1]$ and $5^{2}=25$ quadrature points in each element using the reduced integral (5) and its counter part with $t$ replaced by $-t$. At the focus $x=0$, the results are reported at Table 1, in which the errors are calculated by $e_{k}=|u-u_{GB}|/|u_{GB}|$, and the orders of convergence are obtained by

[TABLE]

We also test the energy errors and orders of convergence at some $t$ in $(0,1)$. The error in energy norm is calculated by $e_{k}={\|u-u_{k}\|_{E}}/{\|u_{k}\|_{E}},$ with $\|v\|_{E}^{2}=\frac{1}{2}\int_{\mathbb{R}^{3}}|u_{t}|^{2}+|\nabla_{x}u|^{2}dx$, evaluated over the ball of radius $r_{1}+t$. The errors and orders of convergence using (5.4) are reported in Table 2.

The numerical results with the gain in the order of convergence agree with the above asymptotic estimate.

The mechanism that leads to a relative error of order $k^{-1}$ in this example is probably the cancelation of terms of order $k^{-1/2}$ in the Gaussian beam superposition in (5). This is the result of the spherical symmetry in this superposition.

6. A 2D example with fold caustics

This section is devoted to the construction of a Gaussian beam superposition with fold caustics for the 2D acoustic wave equation

[TABLE]

Let us consider a gaussian beam superposition with ray paths given by

[TABLE]

where

[TABLE]

These rays are tangent to the unit circle at $(x_{1},x_{2})=(\cos\theta,\sin\theta)$ and propagating in the direction of the tangent $(\sin\theta,-\cos\theta)$. That defines the parameter $\theta$. The parameter $s$ is distance along the ray path, chosen so that $s=0$ on $x_{1}^{2}+x_{2}^{2}=2$ and $s=1$ on $x_{1}^{2}+x_{2}^{2}=1$. More precisely the relation between $r$ and $s$ is (by the Pythagorean Theorem)

[TABLE]

or $s=1-\sqrt{r^{2}-1}$ for $s<1$ and $s=1+\sqrt{r^{2}-1}$ for $s>1$. The phase function associated with these ray paths, which was complicated in euclidian coordinates, is quite simple in $(\theta,s)$. It can be chosen as

[TABLE]

defined for $0\leq s<1$, and $-\pi<\theta<\pi$. The function $\exp(ikS(x_{1},x_{2}))$ will be single-valued on the annulus bounded by the circles of radius $1$ and $\sqrt{2}$ only when $k$ is an integer. In the numerical examples we will take $k$ to be an integer.

The Hessian of $S(x_{1},x_{2})$ has to be a multiple of the orthogonal projection $P^{\perp}$ onto $(\cos\theta,\sin\theta)$, the vector perpendicular to the ray path. We find that at $(x_{1}(0;\theta,s),x_{2}(0;\theta,s))$,

[TABLE]

The Hessian of the phase in the Gaussian beam, $M(t;\theta,s)$, has to be given by

[TABLE]

where $P=I-P^{\perp}$, $\partial_{t}a=0$, $\partial_{t}b+b^{2}=0$, and $(a(0,s),b(0,s))=(i,i+1/(s-1))$. So the Hessian of the phase is a lot like the Hessian in the 3D example. In fact, we have

[TABLE]

The next step in the construction would be to find the amplitude, but for that one needs the phase. That is

[TABLE]

That comes from formulas (1.5) and (1.6) in [9] with one small observation: the function $\phi_{0}(t;z)$ with $z=(\theta,s)$ does not depend on $t$. You can see that from the fact that since the Hamitonaian $H$ is homogeneous of degree one in $p$, $\dot{x}(t;z)\cdot p(t;z)=H(t,x(t;z),p(t;z))$, which forces $\dot{\phi}_{0}(t;z)=0$ (see also equation (3.10c) in [11]).

Continuing, we have $\partial_{t}\phi=-1$ and $\square_{x,t}\phi=-b(t,s)$ when $x=x(t;\theta,s)$. So the solution of the transport equation, $2A_{t}\phi_{t}+(\square_{x,t}\phi)A=0$ is just

[TABLE]

So the complete Gaussian beam superposition will be

[TABLE]

where $1-s$ is the absolute value of the Jacobian of $(x_{1}(0;\theta,s),x_{2}(0,\theta,s))$ with respect to $(\theta,s)$. The contributions from beams built with the other choice, $\partial_{t}\phi=1$, propagate away from the disk $\{|x|\leq 1\}$ as $t$ increases, and are negligible near the caustics on the circle. Hence we have omitted those contributions from all formulas and numerical results below. In the next section we will examine the accuracy of the method numerically.

6.1. Numerical results

In addition to estimates of accuracy in the energy norm, we will also give numerical estimates in the maximum norm. The results in [11] restricted to first order beams with $O(1)$ initial data show that $||u_{GB}(t)-u(t)||_{L^{\infty}}\leq Ck^{-1}$ away from caustics, see [11, estimate (6.1)]. For domains including caustics [11] gives the weaker estimate $||u_{GB}(t)-u(t)||_{L^{\infty}}\leq k^{1/2}$, and these estimates hold in relative norms well. Our numerical results in Table 3 below show that at caustics the order of error is

[TABLE]

with $\alpha(t)$ varying in $(0.5,1)$. We see that the numerical order of error near caustics is greater than the error away from caustics but much smaller than the bound in [11].

We consider the 2D acoustic wave equation (6.1) on $[0,T]\times\Omega$, where $\Omega=[-L/2,L/2]^{2}$ with $L=4$, subject to initial data $(u,\partial_{t}u)(x,0)=(u_{GB},\partial_{t}u_{GB})(x,0)$ and periodic boundary conditions. For the Gaussian beam superposition (6.2) we take initial amplitude

[TABLE]

which is supported on $1+(1-s_{1})^{2}\leq x_{1}^{2}+x_{2}^{2}\leq 1+(1-s_{0})^{2}$ for $s_{0},s_{1}\in(0,1)$, since $1+(s-1)^{2}=r^{2}$.

We use the fast Fourier transform to approximate the “exact solution”, and use it to determine the errors in the Gaussian beam superposition. For $K$ large enough, say $K=1024$, we partition $\Omega$ by a uniform rectangular mesh $\Omega=\left[-L/2:h:L/2-h\right]^{2}$, with $h=L/K$. We obtain the “exact solution” and its derivatives numerically using Matlab 2018a in the following steps.

• •

Step 1 (Initial preparation) We calculate the integrals in the Gaussian beam superposition $u_{GB}(x,t)$; more specifically $u_{GB}(x,0),(u_{GB})_{t}(x,0)$ using integral2 with absolute tolerance $10^{-8}$.

• •

Step 2 (fast Fourier transform) The Fourier transform of (6.1) gives

[TABLE]

where $(\kappa_{1},\kappa_{2})\in\left[\frac{2\pi}{L}(0,1,\cdots,K/2-1,-K/2,-K/2+1,\cdots,-1)\right]^{2}$ is adopted in the fast Fourier transform (fft2) in Matlab.

• •

Step 3 (Solving ODE) The exact solution of (6.3) is determined by

[TABLE]

• •

Step 4 (inverse fast Fourier transform) We obtain the “exact solution” $u$ and its derivatives $\partial_{t}u,\partial_{x_{1}}u,\partial_{x_{2}}u$ through the inverse fast Fourier transform (ifft2 in Matlab) applied to $\widehat{u},\partial_{t}\widehat{u},i\kappa_{1}\widehat{u},i\kappa_{2}\widehat{u}$, respectively.

We test the case $s_{0}=0.25,\ s_{1}=0.75$. With this choice, the wave propagates within the entire computational domain $\Omega$ for $t\leq T=0.8$, and caustics appear only for $0.25=1-s_{1}<t<1-s_{0}=0.75$.

In this example, $u(x,0)=u_{k}(x,0)$. A refined numerical test indicates that

[TABLE]

where the rate $\beta(t)$, shown in Figure 1, is calculated over $N\times N$ meshes with $N=2.5\times 10^{5}$, and of frequencies $k=40960$ and $k=81920$. From this figure, we see that $\beta(t)\sim 0$ when away from caustics, but $\beta(t)$ can go up to about $1/6$ in the presence of caustics.

The experimental orders of convergence (EOC) are obtained by

[TABLE]

where $e_{k}$ is the relative error between the “exact solution” $u(x,t)$ and $u_{k}:=u_{GB}$.

Test case 1. Convergence in $L^{\infty}$ norm.

We first check the $L^{\infty}$ errors and orders of convergence from $t=0.15$ to $t=0.8$. The error in $L^{\infty}$ norm is approximated by

[TABLE]

From the errors and orders of convergence reported in Table 3 obtained using $1024\times 1024$ meshes, we find that the orders of accuracy are decreased in the present of fold caustics. At $t=0.15$ and $0.80$, the orders of accuracy are increased when $k$ increases to $320$, much closer to the desired first order since caustics are not present.

Test case 2. Convergence in energy norm. We next check the energy errors and orders of convergence from $t=0.15$ to $0.8$. The error in energy norm is approximated by

[TABLE]

The results in Table 4 show that $1$st order of accuracy in energy norm is obtained regardless of the appearance of caustics. Note that in contrast to the maximum norm, for energy norm $\|u_{k}(\cdot,t)\|_{E}=\|u(\cdot,0)\|_{E}$.

The gain in order of accuracy in the energy norm indicates the contribution from cancellations of first order beams, this is consistent with the numerical evidence in [9]. However, the order of accuracy in $L^{\infty}$ norm can vary in time due to the presence of caustics; while when away from caustics the uniform first order of accuracy in $L^{\infty}$ norm has been proven in [11].

7. Examples of general superpositions

In this section we discuss the growth rate in $k$ of general superpositions,

[TABLE]

when measured in the energy norm. This will depend on the detailed description of $K_{0}$, and we discuss by examples. To simplify presentation, we only estimate the $L^{2}$-norm of $\nabla_{x}u$ in all examples, instead of computing the whole energy norm. For beams the $L^{2}$-norm of the spatial gradient is always comparable to the $L^{2}$ norm of the initial time derivative.

Let $K_{0}$ be parameterized by $z\in\Sigma$ so that

[TABLE]

Here listed are some typical examples.

**Example 1. ** If the data is concentrated at one point (say, in the case of a point source for stationary problems), one may consider

[TABLE]

In the example presented in §4, we have

[TABLE]

This corresponds to $d=3$ and $m=2$ with

[TABLE]

The asymptotic rate of its energy norm is

[TABLE]

We may also consider the case $m=d$ with

[TABLE]

Let $a(p)$ be a smooth function compactly supported in $p$, and

[TABLE]

Hence $u_{GB}(x,0)=k^{m/2}\hat{a}(kx)e^{-k|x|^{2}/2}$, and

[TABLE]

This implies $\|u_{GB}(\cdot,0)\|_{E}\sim k^{1-\frac{d-m}{4}}$. This together with the result in Theorem 3.1 says that the relative error is no greater than $k^{-N/2}$, as we expected.

**Example 2. ** A more general example of a superposition.

Let $z=(z^{(1)},z^{(2)})$ where $z^{(1)}=(z_{1},...,z_{r})\in\mathbb{R}^{r}$ and $z^{(2)}=(z_{r+1},..,z_{m})\in\mathbb{R}^{m-r}$. Consider the superposition of Gaussian beams in $\mathbb{R}^{d}$

[TABLE]

Here $x^{(3)}=(x_{m+1},...,x_{d})\in\mathbb{R}^{d-m}$. We will take $a(z)=e^{-|z|^{2}/2}$ to make some computations explicit. So $a(z)$ nearly has compact support. We have

[TABLE]

Since

[TABLE]

then

[TABLE]

We have

[TABLE]

This gives

[TABLE]

where $c_{1}$, $c_{2}$ and $c_{3}$ are powers of $2\pi$. The first term in that expression dominates, and we have

[TABLE]

Assuming that the $L^{2}$-norm of $\partial_{t}u_{GB}(x,0)$ is of the same order, we can compare that with $||u(\cdot,t)-u_{GB}(\cdot,t)||_{E}$ for which we have the estimate (for $|t|<T$)

[TABLE]

for first order beams, and get the relative error estimate

[TABLE]

This shows what can happen when initial data is not of form (1.2).

**Example 3. ** For wave equation (1.1) subject to the WKB initial data,

[TABLE]

compactly supported in $\Omega\subset\mathbb{R}^{d}$, one may consider $m=d$ with

[TABLE]

The superposition of the first order Gaussian beam is given by

[TABLE]

where

[TABLE]

with $p_{0}=\nabla_{x}S_{0}(x_{0})$ and $M_{0}=\partial_{x}^{2}S_{0}(x_{0})+iI$. Note that

[TABLE]

Hence the energy norm can be estimated as

[TABLE]

This upper bound is as expected.

**Example 4. ** For the WKB data $e^{(ik-1)|x|^{2}/2}$, we consider

[TABLE]

Note that $|x|^{2}/2=|z|^{2}/2+z\cdot(x-z)+|x-z|^{2}/2$, this superposition corresponds to the case with $p(z)=z$, $x(z)=z$, initial phase $S_{0}(x)=|x|^{2}/2$, and initial amplitude $e^{-|x|^{2}/2}$. Since

[TABLE]

we have

[TABLE]

This implies

[TABLE]

and passing to $\nabla u_{GB}(x,0)$ brings down a factor of order $k$. Hence,

[TABLE]

We may also consider

[TABLE]

where $a$ is assumed to be smooth with compact support. This corresponds to the case with $p(z)=z$, $x(z)=z$, initial phase $S_{0}(x)=|x|^{2}/2$, and initial amplitude $a(x)$, for

[TABLE]

We have

[TABLE]

This implies

[TABLE]

**Example 5. ** This example is a bit surprising.

Let

[TABLE]

In other words $p(z)=z$ and $x(z)=z$. In this case it is easy to compute $u_{GB}(x,0)$. Since

[TABLE]

we have

[TABLE]

We can see this implies

[TABLE]

and passing to $\nabla u_{GB}(x,0)$ brings down factors of $x_{j}$ multiplied by factors of order $k$. Hence,

[TABLE]

Note that here $||u_{GB}(\cdot,0)||_{E}$ is not of order $k$. However, like Example 2, the initial data here is not of form (1.2).

We may consider a more general case in the form

[TABLE]

where $a$ is assumed to be smooth with compact support. Let $\eta=z-x$ and the integral becomes

[TABLE]

Using the Plancherel Theorem one can write

[TABLE]

Now, assuming that $a(z)$ is smooth with compact support, $|\hat{a}(\xi)|\leq C_{N}(1+|\xi|^{2})^{-N}$ for all $N$. So

[TABLE]

Now divide that integral into $I_{1}=\int_{\{|\xi|<k|x|/2\}}$ and $I_{2}=\int_{\{|\xi|>k|x|/2\}}$. Then, taking $N$ large enough that $\int_{\mathbb{R}^{d}}(1+|\xi|^{2})^{-N}d\xi<\infty$, the contribution to $|u_{GB}(x,0)|$ from $I_{1}$ is bounded by

[TABLE]

and the contribution from $I_{2}$ is bounded by

[TABLE]

Finally we split $I_{0}$ into $\chi_{\{|x|>k^{-1/2}\}}(x)I_{0}(x)+\chi_{\{|x|<k^{-1/2}\}}(x)I_{0}(x)$ ($\chi_{E}$ is the characteristic function of $E$) . Using that splitting and taking $N$ sufficiently large ($N=N(M)$), one ends up with for any $M>0$ and $\alpha>1$,

[TABLE]

That leads once more to

[TABLE]

8. Final remarks

We have presented results on superpositions of Gaussian beams of order $N$ in dimension $d$ over arbitrary bounded sets of dimension $m$ in phase space, and shown that the error in the approximation of the exact solution with the same initial data is $O(k^{1-N/2-(d-m)/4})$ in energy norm. This result is sharp for general super-positions. For exact solutions with WKB initial data, i.e. initial data of the form (1.2) our numerical evidence in the case $N=1$ and $d=m$ indicates the stronger estimate $O(1)$, or $O(k^{-1})$ in the relative energy norm as conjectured in [9]. However, the numerical estimates in maximum norm are not uniform in time due to the presence of caustics; while away from caustics we know the relative propagation error in maximum norm is $O(k^{-1})$ as has been proven in [11].

Acknowledgments

This work was supported by the National Science Foundation under Grant RNMS (Ki-Net) 1107291 and by NSF Grant DMS1812666.