Abstract

We consider high-frequency waves satisfying the scalar wave equation with highly oscillatory initial data. The wave speed, and the phase and amplitude of the initial data are assumed to be uncertain, described by a finite number of random variables with known probability distributions. We define quantities of interest (QoIs), or observables, as local averages of the squared modulus of the wave solution. We aim to quantify the regularity of these QoIs in terms of the input random parameters, and the wave length, i.e., to estimate the size of their derivatives. The regularity is important for uncertainty quantification methods based on interpolation in the stochastic space. In particular, the size of the derivatives should be bounded independently of the wave length. In this paper, we are able to show that when these QoIs are approximated by Gaussian beam superpositions, they indeed have this property, despite the highly oscillatory character of the waves.

Keywords

1 Background

Many physical phenomena can be described by propagation of high-frequency waves. Such problems arise, for example, in optics, acoustics, geophysics or oceanography. By high frequency, it is understood that the wavelength is very short compared to the distance travelled by the wave. Uncertainty may often be introduced through one or multiple sources: Either the model parameters are incomplete or it is due to the intrinsic variability of the problem. An example of such uncertain high-frequency wave propagation is an earthquake where seismic waves with uncertain epicenter travel through the layers of the Earth with uncertain soil characteristics.

As a simplified model of the wave propagation, we consider the scalar wave equation

with highly oscillatory initial data represented by the wavelength \(\varepsilon \), which is assumed to be much smaller than both the typical scale of variations in the wave speed and the wave propagation distance. Such initial data generate high-frequency waves propagating in low-frequency media. By \(t\in [0,T]\) we denote the time and by \(\mathbf x = (x_1,\ldots ,x_n)\in \mathbb {R}^n\) the spatial variable. The uncertainty in the model is described by the random variables \(\mathbf{y }= (y_1,\cdots ,y_N) \in \varGamma \subset \mathbb {R}^N\) where N is a positive integer. Their joint probability density is assumed to be bounded. Two sources of uncertainty are considered: the local speed, \(c = c(\mathbf x ,\mathbf y )\), and the initial data, \(B_0=B_0(\mathbf x ,\mathbf y ), B_1=B_1(\mathbf x ,\mathbf y ), \varphi _0=\varphi _0(\mathbf x ,\mathbf y )\). The solution is therefore also a function of the random parameter, \(u^{\varepsilon }=u^{\varepsilon }(t,\mathbf x ,\mathbf y )\). The initial data and the local speed are all assumed to be smooth functions.

In this work, we are interested in observables, or quantities of interest (QoI) defined as functionals of the solution \(u^\varepsilon \) to (1) for each \({\mathbf{y }}\in \varGamma \). These observables often carry important physical information. They are functions of the random variable \(\mathbf{y }\), and we want to determine their statistical moments, for example the expected value and the variance. This means that we first need to solve (1), then evaluate the functional and finally integrate in \(\mathbf y \) over \(\varGamma \), which typically is in a high-dimensional space.

Hence, when computing the statistical moments, we are facing two problems: high frequencies and high-dimensional integrals. Firstly, modeling of a high-frequency problem by direct methods is not feasible. The accuracy of direct methods is determined by the number of points per wavelength. For a given accuracy, the wave resolution cost grows rapidly when shrinking the wavelength. As a consequence, direct numerical methods become very expensive at high frequencies. Secondly, the dimension of the random space is typically large which results in a fast growth of the number of mesh points in standard numerical integration schemes. This phenomenon is referred to as the curse of dimensionality.

In [18], we proposed a method for computing the statistics of quantities of interest that deals efficiently with both difficulties. It consists of two techniques: the Gaussian beam superposition method for treating the high frequencies and the sparse stochastic collocation for treating the uncertainty. The Gaussian beam method is an asymptotic technique largely based on the geometrical optics approach. It has been extensively studied in [4, 16, 23, 29]. The sparse grid collocation method has been used in forward propagation of uncertainty in many PDE models, see [16, 25]. The two methods separately are in use, but their combination offers a new powerful tool.

In order for the sparse stochastic collocation to converge fast, high stochastic regularity of the QoI is crucial. That is, the derivatives of the QoI with respect to \(\mathbf y \) should exist and be small. In the high-frequency setting, the QoI will depend on the wavelength \(\varepsilon \) and it is important that the derivatives are bounded independently of \(\varepsilon \). In particular, the QoI should not oscillate with period \(\varepsilon \). The objective of this paper is to show such bounds for a certain observable.

Natural choices of QoI include the solution itself in a fixed point \((T,\mathbf x _0)\),

In our case, however, neither of these QoIs is interesting due to the high-frequency nature of \(u^{\varepsilon }\). In the first case, the pointwise QoI is itself highly oscillatory in \(\mathbf y \) and its derivatives cannot be bounded by constants that are independent of \(\varepsilon \). In the second case, the mean value tends to zero as \(\varepsilon \rightarrow 0\) due to the stationary phase lemma [9], and little new information is given by the QoI.

Quadratic quantities present another choice of observables. In this paper, we will consider:

where \(\psi \) is again a given real-valued compactly supported smooth function referred to as the QoI test function. QoI (2) has an important physical motivation. Since the square of the wave amplitude \(|u^\varepsilon |^2\) represents energy density, the weighted integration of the wave amplitude represents the weighted average strength of the wave inside the support of \(\psi \). For example, if the solution \(u^{\varepsilon }\) describes the acoustic pressure, \(\mathcal {Q}^\varepsilon \) represents the local acoustic potential energy.

It was conjectured in [18] that the QoI (2) computed by the Gaussian beam method, \({\mathcal {Q}}_{\mathrm{GB}}^\varepsilon \), has high stochastic regularity, more precisely, that

The result (3) may seem counterintuitive as the solution itself oscillates with \(\varepsilon \). However, formal theoretical arguments and numerical examples presented in [18] showed that \(\mathcal {Q}^\varepsilon _{\mathrm{GB}}\) is indeed non-oscillatory. As a consequence, the convergence of the sparse stochastic collocation method was exponential, contrasting the rather slow convergence for Monte Carlo simulations.

The main result of this work is Theorem 1 in Sect. 5, where (3) is proved. This shows the suitability of the sparse grid collocation method for computation of the statistics of the QoI. The proof is based on techniques presented in [16]. We stress that (3) is shown for the approximate solution. Whether it holds also for the exact solution with matching initial data is an open question. A simple argument based on the accuracy of the Gaussian beam solution gives that at least \(|\mathcal {Q}^\varepsilon _{\mathrm{GB}}-\mathcal {Q}^\varepsilon |=O(\varepsilon ^{k/2-1})\), for kth-order beams, see Remark 1. Hence, for sufficiently high-order beams \(\mathcal {Q}^\varepsilon _{\mathrm{GB}}\) is close to \(\mathcal {Q}^\varepsilon \). Based on numerical experiments, we conjecture that it is in fact close for all beam orders.

The paper is organized as follows: In Sect. 2, we introduce precise assumptions. In Sect. 3, we present the sparse collocation technique and a numerical example. Section 4 contains the Gaussian beam approximation description and properties. In Sect. 5, we prove the main result on stochastic regularity of the QoI.

2 Preliminaries and assumptions

In this section, we introduce notation used in the paper and formulate precise assumptions for the high-frequency wave propagation problem (1) with stochastic parameters. We write \(|\mathbf x |\) for the Euclidean norm of a vector \(\mathbf x \in \mathbb {R}^n\). For a multi-index \(\varvec{\alpha }=(\alpha _1,\ldots ,\alpha _n)\in \mathbb {N}^n_0\), we, however, use the standard convention that \(|\varvec{\alpha }|= \alpha _1+\cdots +\alpha _n\). We frequently use the simple estimate,

3 Forward propagation of uncertainty

There are in general two types of numerical methods for propagating uncertainty in PDE models subject to stochastic uncertainty: (1) intrusive methods and (2) non-intrusive methods. Intrusive methods, such as stochastic Galerkin [2, 33], are based on a weak formulation of the problem and are combined with direct simulations such as the finite element method. Consequently, their application to high-frequency waves is not feasible, since direct solvers require a large number of elements to resolve the wave’s high oscillations. On the contrary, non-intrusive methods, such as Monte Carlo [6] and stochastic collocation [1, 20, 32], are sample-based approaches and rely on a set of deterministic models corresponding to a set of realizations. Consequently, asymptotic techniques for approximating deterministic high-frequency waves [5, 22, 27], such as Gaussian beam superposition (presented in Sect. 4), can be employed within the framework of a non-intrusive method.

In this section, we discuss and motivate the applicability and efficiency of two popular non-intrusive methods, Monte Carlo and stochastic collocation, for approximating the statistics of the QoI in (2).

3.1 Non-intrusive numerical methods

Two popular non-intrusive methods, combined with Gaussian beam superposition, have been proposed in the high-frequency wave community: a Monte Carlo method [31] and a sparse stochastic collocation method [18]. In both methods, the solution to problem (1) is approximated in space and time by the Gaussian beam superposition (see Sect. 4). This gives the Gaussian beam QoI \({{\mathcal {Q}}}^\varepsilon _{\text {GB}}(\mathbf{y })\) in (13), which is evaluated at some set of points \(\{ \mathbf{y }^{(m)}\}_{m=1}^{\eta }\). The two methods differ in treating the Gaussian beam QoI in the stochastic space, including the selection of the set of points \(\{ \mathbf{y }^{(m)}\}_{m=1}^{\eta }\). Here, we briefly review and compare the convergence rates of the two methods.

Monte Carlo. The Monte Carlo method approximates the statistical moments of a stochastic quantity by sample averages. For instance, let \(\{ \mathbf{y }^{(m)} \}_{m=1}^\eta \in \varGamma \) be a set of \(\eta \) independent and identically distributed samples of the random vector \(\mathbf{y }\). Then, the Monte Carlo estimator [31] computes the first moment of the quantity \({{\mathcal {Q}}}^\varepsilon \) by sample averages of the Gaussian beam QoI:

for suitable multivariate polynomials, \(\{ L_m (\mathbf{y })\}_{m=1}^{\eta }\), such as Lagrange interpolating polynomials. It is important to note that the statistical moments of the quantity \({{\mathcal {Q}}}^\varepsilon \) can directly be computed without forming the polynomial approximation \( {{\mathcal {Q}}}^\varepsilon _{\text {GB}, \eta } (\mathbf{y })\) in (4). For instance, the first moment of the quantity \({{\mathcal {Q}}}^\varepsilon \) is computed by

The choice of the set of collocation points, \(\{ \mathbf{y }^{(m)} \}_{m=1}^{\eta }\), i.e., the type of computational grid in the N-dimensional stochastic space, is crucial in the computational complexity of stochastic collocation. A full tensor grid, based on the Cartesian product of mono-dimensional grids, can only be used when the dimension of the stochastic space, N, is small, since the computational cost grows exponentially fast with N (the curse of dimensionality). In this case, the grid consists of \(\eta = (p + 1)^N\) collocation points, where p is the polynomial degree in each direction. Alternatively, sparse grids, such as total degree and hyperbolic cross grids, can reduce the curse of dimensionality. In both full tensor and sparse grids, typical choices of grid points to be used include Gauss abscissas [30], Clenshaw–Curtis abscissas [30] and Leja abscissas [24]. The statistical moments of the QoI can then be computed using Gauss or Clenshaw–Curtis quadrature formulas applied to (4). See [18] for more details.

3.2 Convergence of numerical methods

By the central limit theorem [6], if the variance of \(\mathcal {Q}^\varepsilon _{\mathrm{GB}}\) is bounded, then the convergence rate of Monte Carlo estimator is \({{\mathcal {O}}}(\eta ^{-1/2})\). Therefore, while being very flexible and easy to implement, Monte Carlo sampling features a very slow convergence rate. Multi-level Monte Carlo [7, 19], multi-index Monte Carlo [8], and quasi-Monte Carlo methods [11] have been proposed to accelerate the convergence of Monte Carlo sampling. It is to be noted that these recent techniques have not been applied to high-frequency wave propagation problems.

The rate of convergence of stochastic collocation for approximating a stochastic function \(v(\mathbf{y })\) depends on the regularity of the function in stochastic space, i.e., the regularity of the mapping \(v: \varGamma \rightarrow {{\mathbb {R}}}\). Fast convergence with respect to the number of collocation points is attained in the presence of high stochastic regularity. For instance, for a stochastic function \(v(\mathbf{y })\) with \(s \ge 1\) bounded mixed derivatives, \(\partial _{y_1}^{s} \, \partial _{y_2}^{s} \cdots \partial _{y_N}^{s} v \in L^\infty (\varGamma )\), the error in stochastic collocation on sparse Smolyak grids based on Gauss–Legendre abscissas satisfies (see Theorem 5 in [20])

Here, the constant C depends on \(s, N, \rho \) and the function v, but not on \(\eta \). Therefore, the larger s, or equivalently the higher the regularity of v, the faster the convergence rate in \(\eta \). In particular, when the stochastic function has infinitely many bounded mixed derivatives, the convergence rate is close to spectral. See also [3, 20, 21].

In the case of high-frequency waves with a highly oscillatory wave solution, the constant C depends in general also on the wavelength \(\varepsilon \). To maintain fast convergence, it is therefore crucial that C remains bounded as \(\varepsilon \rightarrow 0\). Since C depends on the size of the \(\mathbf y \)-derivatives of v on \(\varGamma \), we need uniform bounds for the \(\mathbf{y }\)-derivatives of \(\mathcal {Q}^\varepsilon _{\text {GB}}(\mathbf y )\) independent of \(\varepsilon \) to obtain fast convergence for all \(\varepsilon \). We emphasize that since the Gaussian beam quantity \(\mathcal {Q}^\varepsilon _{\text {GB}}(\mathbf y )\) is used in the stochastic collocation approximation (4)–(5), it is the stochastic regularity of the Gaussian beam quantity \(\mathcal {Q}^\varepsilon _{\text {GB}}(\mathbf y )\) that determines the convergence rate, not the regularity of the analytical quantity \(\mathcal {Q}^\varepsilon (\mathbf y )\).

3.3 Example

Let us present a numerical example to illustrate the discussion in Sect. 3.2. We introduce two 2D pulses, nearly disjoint in the beginning and moving toward each other. At a finite time, they start to overlap which leads to oscillations of order \(\varepsilon \) in the magnitude of the solution. We consider three uniformly distributed independent random variables \(\mathbf{y }=(y_1, y_2, y_3)\) characterizing a random speed, \(c=c(\mathbf y )=y_1\), and a random initial position of the second pulse, \(\mathbf s _2 = \mathbf s _2(\mathbf y ) = (y_2,y_3)\). We use the smooth test function

The amplitude \(B_1\) is chosen such that the pulses move toward each other, without producing any backward propagating waves. The first-order Gaussian beam method from Sect. 4 has been utilized to compute the solution \(u^{\varepsilon }\). Figure 1 shows the absolute value of \(u^{\varepsilon }\) with \(\mathbf y =(0.8,1,0.5)\) and wavelength \(\varepsilon = 1/40\) at various times. At time \(T = 1\), the pulses reach the central circle indicating the support of the QoI test function with oscillations emerging in the overlap area.

Fig. 1

Time evolution of \(|u^{\varepsilon }(T,\mathbf{x },\mathbf{y })|\) with wavelength \(\varepsilon = 1/40\) and \(\mathbf y =(0.8,1,0.5)\) in the example in Sect. 3.3. Panels show snapshots of \(|u^\varepsilon |\) at the given time values T. The central circle indicates the support of the QoI test function

In Fig. 2, we plot the quantity of interest (13) and its first and second derivatives along the line \(\mathbf{y }(r) = (0.8+0.2r,1+0.5r,0.5r), r\in [0,1]\) at time \(T=1\). We note that the solution \(u^{\varepsilon }\) oscillates with \(\varepsilon \); however, this is not the case for the QoI and its derivatives. This suggests bounds of the type (3).

versus the number of collocation points \(\eta \) for various wavelengths. The Gaussian beam quantity is computed by (4) on the sparse Smolyak grid at time \(T=1\). To obtain a reference solution, we use \(\eta _{\text {ref}}\) = 30,465 collocation points. The figure also compares convergence of the Monte Carlo and the stochastic collocation methods. A simple linear regression through the Monte Carlo data points shows that the rate of convergence of Monte Carlo is approximately 0.46. We observe a fast spectral convergence of the stochastic collocation error due to the high stochastic regularity of the Gaussian beam QoI. Moreover, the decay rate does not deteriorate as \(\varepsilon \) decreases.

4 Gaussian beams

In this section, we briefly describe the Gaussian beam approximation. We restrict the description to the points that are relevant for the analysis in Sect. 5. For a more detailed account with a general derivation for hyperbolic equations, dispersive wave equations and Helmholtz equation, we refer to [12–14, 16, 26, 28].

Individual Gaussian beams are asymptotic solutions to the wave equation that concentrate around a central ray in space-time. We denote the k-th-order Gaussian beam and the central ray starting at \(\mathbf z \in K_0\) by \(v_k(t,\mathbf y ;\mathbf z )\) and \(\mathbf q (t,\mathbf y ;\mathbf z )\), respectively. The beam has the following form,

where the integration in \(\mathbf z \) is over the support of the initial data \(K_0\subset \mathbb {R}^n\). Since the wave equation is linear, the superposition of beams is still an asymptotic solution. The function \(\varrho _\eta \in C^\infty (\mathbb {R}^n)\) is a real-valued cutoff function with radius \(0<\eta \le \infty \) satisfying,

As shown in (P4), if \(\eta >0\) is sufficiently small, it is ensured that \(\text {Im}\ \varPhi _k>0\) on the support of \(\varrho _\eta \) and the Gaussian beam superposition is well behaved. For first-order beams, \(k=1\), the cutoff function is not needed and we can take \(\eta =\infty \).

Since the wave equation (1) is a second-order equation, two modes and two Gaussian beam superpositions are needed, one for forward and one for backward propagating waves. We denote the corresponding coefficients by \({}^+\) and \({}^-\), respectively, and write the solution

Each Gaussian beam \(v^\pm _{k}(t,\mathbf x ,\mathbf y ;\mathbf z )\) requires initial values for the central ray and all of the amplitude and phase Taylor coefficients. The appropriate choice of these initial values will make \({\tilde{u}}_k(0,\mathbf x ,\mathbf y )\) asymptotically converge to the initial conditions in (1). As shown in [16], initial data for the central ray and phase coefficients should be chosen as follows:

The initial data for the amplitude coefficients are more complicated, and we refer the reader to [16] for more detailed description.

In the remainder of the paper, we will restrict ourselves to only one of the two modes and drop the ± subscript. This means that we only consider high-frequency wave solutions that propagate in one direction from the initial source. Note that the interaction between the two modes could potentially introduce more oscillations in \(\mathcal {Q}^\varepsilon _{\mathrm{GB}}\). An analysis of this will be reported elsewhere.

4.1 Gaussian beam quantity of interest

Instead of the full quantity of interest defined in (2), we here consider the corresponding quantity of interest computed based on the Gaussian beam approximation with one mode,

where \(u_{\mathrm{k}}\) is given by (9) and \(v_k\) is either \(v_k^+\) or \(v_k^-\), with corresponding \(\mathbf q ^+\) or \(\mathbf q ^-\) respectively. Since \(u_k\) is a good approximation of \(u^\varepsilon \), the QoI \(\mathcal Q^\varepsilon _{GB}\) in (13) is a good approximation to \(\mathcal Q^\varepsilon \), see Remark 1.

if initial data are sufficiently well approximated by the beams. This immediately gives the same estimate for \({{\mathcal {Q}}}_{\mathrm{GB}}^\varepsilon \) when \(\psi \ge 0\). Indeed, for a fixed \({\mathbf{y }}\),

Below we show that \({{\mathcal {Q}}}^\varepsilon _{\mathrm{GB}}\) is bounded independently of \(\varepsilon \), and the estimate therefore implies that \({{\mathcal {Q}}}^\varepsilon _{\mathrm{GB}}({\mathbf{y }})\) is close to \({{\mathcal {Q}}}^\varepsilon ({\mathbf{y }})\) for small \(\varepsilon \), at least if we use beams of third order or higher. We do not believe this estimate is sharp, however. In particular, in computations we see a small difference already with first-order beams.

4.2 Gaussian beam phase properties

In this section, we present some consequences of assumptions (A1)–(A4) for the Gaussian beam coefficient functions, the whole phase \(\varPhi _k\) and the symmetrized phase \(\varPsi \). The properties are of great importance in the subsequent estimates in Sect. 5. They also show the significance of the cutoff function introduced in Sect. 4 and indicate how the width \(\eta \) should be chosen. Following (P4), we make the definition:

Definition 1

The cutoff width \(\eta \) used for the Gaussian beam approximation is called admissible for \(\varPhi _k\) if it is small enough in the sense of property (P4) in Proposition 1.

Throughout this paper, we will consider the cutoff function \(\varrho _\eta \) in the Gaussian beam approximation (9) to be chosen such that \(\eta \) is admissible in the sense of Definition 1.

Because of the cutoff function, the terms \(g_{j,m,{\varvec{\alpha }},{\varvec{\beta }}}\) in (18) are zero unless \(|\mathbf x -\varDelta \mathbf q |\) and \(|\mathbf x +\varDelta \mathbf q |\) are both small enough. We introduce the subsets \(\varOmega _\mu \subset \mathbb {R}^n\), parameterized by \(\mu \), to describe this

The proof of these properties essentially follows the same steps as in the corresponding lemmas in [15, 16]. The details are provided in [17].

5 Stochastic regularity

In this section, we prove the main result of the paper, namely that the derivative \(\partial _{\mathbf{y }}^{\varvec{\varvec{\sigma }}} {\mathcal {Q}}_{\mathrm{GB}}^\varepsilon \) is uniformly bounded in \(\mathbf{y }\) and \(\varepsilon \). More precisely, we show the following theorem.

Theorem 1

Under assumptions (A1)–(A6), for the Gaussian beam quantity of interest (13) where the cutoff parameter \(\eta \) is admissible in the sense of Definition 1, there are constants \(C_\sigma \), independent of \(\varepsilon \), such that

In Sect. 5.1, we prove various prerequisite lemmas needed for the proof of Theorem 1 and introduce notation used in the rest of the paper. Section 5.2 is devoted to rewriting the QoI and its derivatives. In Sect. 5.3, several important estimates are presented. In Sect. 5.4, the actual proof of Theorem 1 is carried out.

In this section, we first rewrite the integral \(I^\varepsilon \) (Proposition 2) employing Lemma 3, and phase \(\varPsi \) (Lemma 4). The form of \(\partial ^{\varvec{\sigma }}_\mathbf{y }I^\varepsilon \) then follows from Proposition 3.

where \(\bar{g}_{j,\varvec{\alpha },\varvec{\beta }}\in \mathcal U_\eta \) are either sums of \(\bar{g}_{j,m,\varvec{\alpha },\varvec{\beta }}\in \mathcal U_\eta \) in (18) or identically zero. Recalling Lemma 3, this can be reformulated as

where \(m_{\varvec{\alpha }}, \ell _{\varvec{\alpha },\varvec{\beta }} \in \mathcal V\). Again, we used (P1) and (24). The term \((\mathbf x + \varDelta \mathbf q )^T{M^*(\mathbf y ;\mathbf z )}(\mathbf x +\varDelta \mathbf q )\) can be formulated in a similar fashion. The last two terms in (31) include the quantity \((\mathbf x \pm \varDelta \mathbf q )^{\varvec{\gamma }}\), which can be expanded by (27), namely

for all \(\varvec{\varvec{\sigma }} \in {\mathbb {N}}^n_0\). First, we show that the statement (32) is valid for \(\varvec{\varvec{\sigma }} \equiv 0\). Next, while assuming that the claim holds for \(\varvec{\varvec{\sigma }} \in {\mathbb {N}}^n_0\) we prove that it is valid for any \(\bar{\varvec{\varvec{\sigma }}} = \varvec{\varvec{\sigma }} + \mathbf{e }_\ell \).

1.

\(\varvec{\varvec{\sigma }}= 0\): The sums in (29) are (32) in this case, with \(b_{j,\varvec{\alpha },\varvec{\beta },0}=g_{j,\varvec{\alpha },\varvec{\beta }}\).

where we recalled (35) in the last inequality. We note that the constant C may depend on n and \(\varvec{\beta }\), but is independent of \(\mathbf z , \mathbf z ', \mathbf{y }\) and \(\varepsilon \). \(\square \)

Proof

Take \(\mu \) and \(\theta \) such that they are small enough in the sense of (P6), and also \(\mu <\eta \). We simplify our notation by omitting the functions’ arguments and introduce two smooth test functions \(\varphi _1, \varphi _2 \in C_{\text {c}}^{\infty }({{\mathbb {R}}}^n)\) with the following properties:

We note that the constant C may depend on \((n, \varvec{\beta })\), but is independent of \(\mathbf{y }\) and \(\varepsilon \). This proves the needed estimate for \({I}_1\) for the case \(K=0\) as well as the case \(K>0\) when \(\mathbf z =\mathbf z '\). Therefore, in the remainder of the proof we will consider the case \(K \ne 0\) and \(\mathbf z \ne \mathbf z '\).

From (P6), since \(\mu \) and \(\theta \) are small enough, there exists a constant \(C>0\) such that \(\inf _\mathbf{x \in \varOmega _{\mu }} |\nabla \varPsi |\ge C|\mathbf z -\mathbf z '|>0\). Since the gradient of \(\varPsi \) does not vanish, we employ the non-stationary phase lemma (see [10]). Since \(g_1\in \mathcal U_\mu \), we get

where the last inequality follows because \(g_1\) has bounded derivatives with respect to \(\mathbf{x }\), independent of \(\mathbf{y }\) and \(\varepsilon \), which is a consequence of (A4) and the smoothness of \(\varphi _1\) and g on the compact sets \(\varOmega _{\mu }\) and \(\varOmega _{\eta }\), respectively. Using (P4) and (35), we have

Declarations

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Authors’ Affiliations

(1)

Department of Mathematics and Swedish e-Science Research Center (SeRC), KTH

(2)

Department of Mathematics and Statistics, The University of New Mexico