A fast direct imaging method for the inverse obstacle scattering problem with nonlinear point scatterers

TL;DR

This paper introduces a fast, FFT-based direct imaging method for inverse obstacle scattering that leverages nonlinear point scatterers to generate high harmonics, improving imaging resolution of extended obstacles.

Contribution

A novel imaging approach that incorporates nonlinear point scatterers and high harmonic generation within a generalized Foldy-Lax framework for better obstacle reconstruction.

Findings

Effective imaging of extended obstacles demonstrated through numerical experiments.

Enhanced resolution achieved using nonlinear point scatterers and high harmonic generation.

Abstract

Consider the scattering of a time-harmonic plane wave by heterogeneous media consisting of linear or nonlinear point scatterers and extended obstacles. A generalized Foldy-Lax formulation is developed to take fully into account of the multiple scattering by the complex media. A new imaging function is proposed and an FFT-based direct imaging method is developed for the inverse obstacle scattering problem, which is to reconstruct the shape of the extended obstacles. The novel idea is to utilize the nonlinear point scatterers to excite high harmonic generation so that enhanced imaging resolution can be achieved. Numerical experiments are presented to demonstrate the effectiveness of the proposed method.

Tables

Table 1. Table 1. Parameters used in the numerical experiments.

$N_{\mathrm{point}}$	number of point scatterers
$N_{\mathrm{boundary}}$	number of points to discretize the boundary of the extended scatterer(s)
$N_{\mathrm{direction}}$	number of incident and observation directions
$N_{\mathrm{sampling}}$	number of sampling points along the $x$- and $y$-direction
$T_{\mathrm{invert}}$	time (in seconds) to invert (factorize) the scattering matrix
$T_{\mathrm{solver}}$	time (in seconds) to solve the linear system for one incidence
$T_{\mathrm{ffp}}$	time (in seconds) to evaluate the far-field patterns
$T_{\mathrm{NUFFT}}$	time (in seconds) to apply the NUFFT to evaluate the imaging function

Table 2. Table 2. Time (in seconds) to solve the linear system ( 6.2 ) on an HP workstation.

$N_{\mathrm{point}}$	$N_{\mathrm{boundary}}$	Method 1	Method 2	Method 3
$1000$	$600$	$0.16$	$1.42$	$0.89$
$10000$	$600$	$8.9$	fail to converge	fail to converge

Table 3. Table 3. Results for imaging two extended scatterers surrounded by linear point scatterers.

	$\kappa$	$N_{\mathrm{point}}$	$N_{\mathrm{boundary}}$	$N_{\mathrm{direction}}$	$N_{\mathrm{sampling}}$	$T_{\mathrm{invert}}$	$T_{\mathrm{sampling}}$	$T_{\mathrm{ffp}}$	$T_{\mathrm{NUFFT}}$
Example 1	10	1000	600	360	500	7.95e-2	2.56e-3	2.23e-2	2.46e-1
Example 2	50	1000	4800	1800	500	1.61	2.17e-2	3.49e-1	3.70

Table 4. Table 4. Results for imaging the extended scatterers surrounded by quadratically nonlinear point scatterers.

	$\kappa$	$N_{\mathrm{point}}$	$N_{\mathrm{boundary}}$	$N_{\mathrm{direction}}$	$N_{\mathrm{sampling}}$	$T_{\mathrm{invert}}$	$T_{\mathrm{solver}}$	$T_{\mathrm{ffp}}$	$T_{\mathrm{NUFFT}}$
Example 3	2	2	600	360	500	1.22e-3	8.52e-3	1.39e-2	3.39e-1
Example 4	5	2	1200	360	500	1.90e-1	3.99e-3	1.96e-2	3.84e-1

Table 5. Table 5. Results for imaging the extended scatterers surrounded by cubically nonlinear point scatterers.

	$\kappa$	$N_{\mathrm{point}}$	$N_{\mathrm{boundary}}$	$N_{\mathrm{direction}}$	$N_{\mathrm{sampling}}$	$T_{\mathrm{invert}}$	$T_{\mathrm{solver}}$	$T_{\mathrm{ffp}}$	$T_{\mathrm{NUFFT}}$
Example 5	2	2	600	360	500	1.24e-3	1.00e-2	1.22e-2	4.33e-1
Example 6	5	2	1200	360	500	4.22e-3	1.91e-2	2.11e-2	4.41e-1

Equations

$${\varphi }_{\mathrm{inc}}(\mathbf{r})=e^{\mathrm{i}\kappa \mathbf{r}\cdot \mathbf{d}},\mathbf{r}\in {\mathbb{R}}^2,$$

$$\Delta {\varphi }_{\mathrm{inc}}+{\kappa }^2{\varphi }_{\mathrm{inc}}=0\text{in}{\mathbb{R}}^2.$$

$$\Delta \varphi (\mathbf{r})+{\kappa }^2\varphi (\mathbf{r})=-\sum _{k=1}^m{\sigma }_k{\varphi }_k\delta (\mathbf{r}-{\mathbf{r}}_k),\mathbf{r}\in {\mathbb{R}}^2,$$

$$\varphi ={\varphi }_{\mathrm{inc}}+\psi .$$

$$\Delta \psi (\mathbf{r})+{\kappa }^2\psi (\mathbf{r})=-\sum _{k=1}^m{\sigma }_k{\varphi }_k\delta (\mathbf{r}-{\mathbf{r}}_k),\mathbf{r}\in {\mathbb{R}}^2.$$

$$\underset{r\to \infty }{\lim }{r}^{1/2}({\partial }_r\psi -\mathrm{i}\kappa \psi )=0,r=|\mathbf{r}|.$$

$$\psi (\mathbf{r})=\sum _{k=1}^m{\sigma }_k{\varphi }_k{G}_{\kappa }(\mathbf{r},{\mathbf{r}}_k),$$

$$G_{\kappa }(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{\mathrm{i}}{4}{H}_0^{(1)}(\kappa |\mathbf{r}-{\mathbf{r}}^{\prime }|)$$

$$\varphi (\mathbf{r})={\varphi }_{\mathrm{inc}}(\mathbf{r})+\sum _{k=1}^m{\sigma }_k{\varphi }_k{G}_{\kappa }(\mathbf{r},{\mathbf{r}}_k).$$

$${\varphi }_i-\sum _{\stackrel{k=1}{k\ne i}}^m{\sigma }_k{\varphi }_k{G}_{\kappa }({\mathbf{r}}_i,{\mathbf{r}}_k)={\varphi }_{\mathrm{inc}}({\mathbf{r}}_i),$$

$$\Delta {\varphi }^{(1)}(\mathbf{r})+{\kappa }_1^2{\varphi }^{(1)}(\mathbf{r})$$

$$\Delta {\varphi }^{(2)}(\mathbf{r})+{\kappa }_2^2{\varphi }^{(2)}(\mathbf{r})$$

$${\varphi }^{(1)}={\varphi }_{\mathrm{inc}}+{\psi }^{(1)},{\varphi }^{(2)}={\psi }^{(2)},$$

$$\Delta {\psi }^{(1)}(\mathbf{r})+{\kappa }_1^2{\psi }^{(1)}(\mathbf{r})$$

$$\Delta {\psi }^{(2)}(\mathbf{r})+{\kappa }_2^2{\psi }^{(2)}(\mathbf{r})$$

$$\underset{r\to \infty }{\lim }{r}^{1/2}\left({\partial }_r{\psi }^{(j)}-\mathrm{i}{\kappa }_j{\psi }^{(j)}\right)=0,r=|\mathbf{r}|.$$

$${\psi }^{(1)}(\mathbf{r})$$

$${\psi }^{(2)}(\mathbf{r})$$

$${\varphi }^{(1)}(\mathbf{r})$$

$${\varphi }^{(2)}(\mathbf{r})$$

$${\varphi }_i^{(1)}$$

$${\varphi }_i^{(2)}$$

$$\Delta {\varphi }^{(1)}(\mathbf{r})+{\kappa }_1^2{\varphi }^{(1)}(\mathbf{r})$$

$$\Delta {\varphi }^{(3)}(\mathbf{r})+{\kappa }_2^2{\varphi }^{(3)}(\mathbf{r})$$

$${\varphi }^{(1)}={\varphi }_{\mathrm{inc}}+{\psi }^{(1)},{\varphi }^{(3)}={\psi }^{(3)},$$

$$\Delta {\psi }^{(1)}(\mathbf{r})+{\kappa }_1^2{\psi }^{(1)}(\mathbf{r})$$

$$\Delta {\psi }^{(3)}(\mathbf{r})+{\kappa }_2^2{\psi }^{(3)}(\mathbf{r})$$

$$\underset{r\to \infty }{\lim }{r}^{1/2}\left({\partial }_r{\psi }^{(j)}-\mathrm{i}{\kappa }_j{\psi }^{(j)}\right)=0,r=|\mathbf{r}|.$$

$${\psi }^{(1)}(\mathbf{r})$$

$${\psi }^{(3)}(\mathbf{r})$$

$${\varphi }^{(1)}(\mathbf{r})$$

$${\varphi }^{(3)}(\mathbf{r})$$

$${\varphi }_i^{(1)}$$

$${\varphi }_i^{(3)}$$

$$\Delta \varphi +{\kappa }^2\varphi =0\text{in}{\mathbb{R}}^2\setminus \bar{D}.$$

Full text

The work of J. Lai was supported in part by the Funds for Creative Research Groups of NSFC (No. 11621101) and the Major Research Plan of NSFC (No. 91630309). The research of M. Li was supported partially by the National Youth Science Foundation of China (Grant no. 11401423). The research of P. Li was supported in part by the NSF grant DMS-1151308.

1. Introduction

In scattering theory, one of the basic problems is the scattering of a time-harmonic plane wave by an impenetrable medium, which is called the obstacle scattering problem [19]. Given the incident wave, the direct obstacle scattering problem is to determine the scattered wave for the known obstacle; while the inverse obstacle scattering problem is to reconstruct the shape of the obstacle from either near-field or far-field measurement of the scattered wave. The obstacle scattering problem has played a fundamental role in diverse scientific areas such as geophysical exploration, radar and sonar, medical imaging, and nondestructive testing.

The inverse obstacle scattering problem is challenging due to nonlinearity and ill-posedness. It has been extensively studied by many researchers. Various computational methods have been developed to overcome the issues and solve the inverse problem. Broadly speaking, these methods can be classified into two types: optimization based iterative methods [5] and imaging based direct methods [1, 2, 3, 12, 16, 17, 18, 21, 22, 27, 28, 29, 30, 31, 35, 39, 40, 43]. The former are known as quantitative methods and aim at recovering the functions which parameterize the obstacles. The latter are usually called qualitative methods and attempt to design imaging functions which highlight the obstacles. According to the Rayleigh criterion, there is a resolution limit, roughly one half the wavelength, on the accuracy of the reconstruction for a given incident wave [8, 9, 44]. To improve the resolution, one approach is simply to use an incident wave with shorter wavelength or higher frequency as an illumination. A topical review can be found in [7] on computational approaches and mathematical analysis for solving inverse scattering problems by multi-frequencies.

In this paper, we consider the inverse obstacle scattering problem with an emphasis on resolution enhancement. Motivated by nonlinear optics [10, 11, 38], we utilize nonlinear point scatterers to excite high harmonic generation so that enhanced resolution can be achieved to reconstruct the obstacles. To realize this idea, we have to consider the scattering problem of a time-harmonic plane incident wave by a heterogeneous medium, which consists of small scale point scatterers and wavelength comparable extended obstacles. The main purpose of this work is threefold:

1. (1)

develop a generalized Foldy–Lax formulation and design an efficient solver for the scattering problem involving point scatterers and extended obstacles; 2. (2)

propose a new imaging function and develop an FFT-based method to efficiently evaluate the imaging function on large sampling points; 3. (3)

explore the features of high harmonic generation for nonlinear optics and apply them to the area of inverse scattering to achieve enhanced imaging resolution.

The Foldy–Lax formulation was developed in [23, 41] to describe the scattering of an incoming wave by a group of linear point scatterers. The scattered field can be computed by solving a self-contained system of liner equations. Using a unified approach, we first extend the Foldy–Lax formulation to handle a group of linear, quadratically nonlinear, or cubically nonlinear point scatterers. It is known that the Foldy–Lax formulation is only appropriate for media whose sizes are much smaller than the wavelength [13, 14, 15, 20, 42]. It is no longer adequate for the scattering by wavelength comparable media [24, 37]. The boundary integral equation method is particularly useful for the scattering by extended obstacles. The scattering problem becomes sophisticated when both the point scatterers and extended obstacles are present, as seen in Figure 1. The generalized Foldy–Lax formulation has been studied in [6, 32, 33, 34] to take into account the multiple scattering between linear point scatterers and extended obstacles. Here we develop a generalized Foldy–Lax formulation to fully take into account the multiple scattering between the nonlinear point scatterers and extended obstacles. The generalized formulation couples the Foldy–Lax and boundary integral equation formulations and is self-contained linear or nonlinear system of equations. The coupled system needs to be solved numerically in order to obtain the scattered field and its far-field pattern. For linear point scatterers, we apply LU factorization based direct solver to solve the coupled linear system of equations; for nonlinear point scatterers, we propose an efficient nonlinear solver which combines the Schur complement technique and trust region Newton type method.

The imaging based methods do not require solving any direct problem, which make them very attractive to solve the inverse scattering problem. But, they could be still very time-consuming when evaluating the imaging functions on large sampling points. The reason is that the evaluation usually involves the dense matrix-vector multiplication at each sampling point. To overcome this issue, we construct two illumination vectors and propose a new imaging function. The function, under an appropriate transformation, can be taken as the Fourier transform of the response matrix from the frequency space into the physical space. However, the frequency sampling points are not uniform, which implies the standard FFT can not be applied. We propose an acceleration technique based on the non-uniform fast Fourier transform (NUFFT)[25], which highly reduces the computational complexity and accelerates the evaluation.

Using the quadratically or cubically nonlinear point scatterers, the second or third harmonic generation is excited to image the extended obstacles. Essentially, this approach is equivalent to using double or triple frequency wave to illuminate the obstacles. As a consequence, enhanced resolution can be achieved. However, we find out that the location of the nonlinear point scatterers is crucial to the imaging. Interestingly, they should be aligned up in the same direction as the propagation direction of the plane incident wave so that the correct reconstruction results can appear with the desirable imaging resolution. Numerical experiments are shown to demonstrate the effectiveness of the proposed method. To make the paper self-contained, we briefly introduce in the appendix the nonlinear wave equations and nonlinear point scatterer models.

The paper is organized as follows. In section 2, the Foldy–Lax formulation is introduced for the scattering by a group of point scatterers. In section 3, the boundary integral formulation is briefly reviewed for the scattering by extended scatterers. Section 4 introduces the Generalized Foldy–Lax formulation for the scattering by mixed scatterers. Section 5 is devoted to the inverse scattering problem, where the fast direct imaging method is developed. Numerical experiments are given in section 6. The paper is concluded in section 7.

2. Foldy–Lax formulation

In this section, we introduce the Foldy–Lax formulation for the scattering of a plane incident wave by a group of linear, quadratically nonlinear, and cubically nonlinear point scatterers.

2.1. Linear point scatterers

Consider a collection of $m$ separated linear and isotropic point scatterers located at ${\bf r}_{k}\in\mathbb{R}^{2},k=1,\dots,m$. Let $\phi_{\rm inc}$ be the plane incident wave given explicitly by

[TABLE]

where $\kappa=\omega/c$ is the wavenumber, $\omega>0$ is the angular frequency and $c$ is the speed of wave propagation, and ${\bf d}\in\mathbb{S}^{1}$ is the unit propagation direction. It is easy to verify the incident field satisfies the Helmholtz equation:

[TABLE]

As is shown in (B.2), the total field $\phi$ satisfies

[TABLE]

where $\sigma_{k}>0$ is the scattering coefficient for the $k$-th point scatterer, $\phi_{k}$ is the external field acting on the $k$-th point scatterer, and $\delta$ is the Dirac delta function. We comment that $\phi_{k}$ is also the total field at the location ${\bf r}_{k}$, i.e., $\phi_{k}=\phi({\bf r}_{k})$.

The total field $\phi$ consists of the incident field $\phi_{\rm inc}$ and the scattered field $\psi$:

[TABLE]

Subtracting (2.2) from (2.3) yields

[TABLE]

The scattered field is required to satisfy the Sommerfeld radiation condition:

[TABLE]

It follows from (2.4) and (2.5) that the scattered field can be written as

[TABLE]

where

[TABLE]

is the free space Green’s function for the two-dimensional Helmholtz equation. Here $H_{0}^{(1)}$ is the Hankel function of the first kind with order zero. It is left to compute the external fields $\phi_{k}$ in order to compute the scattered field (2.6).

Adding the incident field on both sides of (2.6) gives

[TABLE]

Evaluating (2.7) at ${\bf r}_{i}$ and excluding the self-interaction to avoid the singularity of the Green function, we obtain a linear system of algebraic equations for $\phi_{k}$:

[TABLE]

which is known as the Foldy–Lax formulation.

2.2. Quadratically nonlinear point scatterers

We assume that the nonlinear point scatterers respond to the incoming wave quadratically and the nonlinearity is weak. Let the plane incident wave $\phi_{\rm inc}$, given in (2.1), be of a single frequency $\omega$. The point scatterers generate fields at frequencies $\omega_{1}=\omega$ and $\omega_{2}=2\omega$ due to the quadratic nonlinearity, which is known as the second harmonic generation. Let $\kappa_{j}=\omega_{j}/c$ and $\phi^{(j)}$ be the field at the frequency $\omega_{j},j=1,2$.

As is known in (B.3), the fields $\phi^{(j)}$ satisfy the coupled Helmholtz equations in $\mathbb{R}^{2}$:

[TABLE]

where the bar denotes the complex conjugate, $\sigma_{k,1}^{(1)}$ and $\sigma_{k,2}^{(1)}$ are the linear scattering coefficients for the $k$-th point scatterer, $\sigma_{k,1}^{(2)}$ and $\sigma_{k,2}^{(2)}$ are the quadratically nonlinear scattering coefficients for the $k$-th point scatterer, and $\phi_{k}^{(j)}$ is the external field acting on the $k$-th point scatterer at frequency $\omega_{j}$. The fields $\phi^{(j)}$ satisfy the following relationship:

[TABLE]

where $\psi^{(j)}$ is the scattered field corresponding to the wavenumber $\kappa_{j}$. It can be verified that the scattered fields satisfy

[TABLE]

In addition, they are required to satisfy the Sommerfeld radiation condition

[TABLE]

It follows from (2.10) that the scattered fields satisfy

[TABLE]

where $G_{\kappa_{j}}$ is the free space Green’s function for the two-dimensional Helmholtz equation at the wavenumber $\kappa_{j}$. Adding the incident field on both sides of (2.11a) and noting $\phi^{(2)}=\psi^{(2)}$ for (2.11b), we obtain

[TABLE]

Similarly, evaluating (2.12) at ${\bf r}_{i}$ and excluding the self-interaction yield a nonlinear system of equations for $\phi_{k}^{(j)}$:

[TABLE]

which is the Foldy–Lax formulation for point scatterers with quadratic nonlinearity.

2.3. Cubically nonlinear point scatterers

Taking the plane incident wave $\phi_{\rm inc}$ with frequency $\omega$, we consider the scattering by point scatterers with weak cubic nonlinearity. The interaction gives rise to fields with frequencies $\omega_{1}=\omega$ and $\omega_{3}=3\omega$, which is called the third harmonic generation. Let $\kappa_{j}=\omega_{j}/c$ and denote the field at frequency $\omega_{j}$ by $\phi^{(j)},j=1,3$.

By (B.4), the fields $\phi^{(j)}$ satisfy the following coupled Helmholtz equations in $\mathbb{R}^{2}$:

[TABLE]

where $\sigma_{k,1}^{(1)}$ and $\sigma_{k,2}^{(1)}$ are the linear scattering coefficients for the $k$-th point scatterer, $\sigma_{k,1}^{(3)}$, $\sigma_{k,2}^{(3)}$ and $\sigma_{k,3}^{(3)}$ are the cubically nonlinear scattering coefficients for the $k$-th point scatterer, and $\phi_{k}^{(j)}$ is the external field acting on the $k$-th point scatterer at frequency $\omega_{j}$. The fields $\phi^{(j)}$ satisfy the following relationship:

[TABLE]

where $\psi^{(j)}$ is the scattered field corresponding to the wavenumber $\kappa_{j}$. It can be verified that the scattered fields satisfy

[TABLE]

In addition, they are required to satisfy the Sommerfeld radiation condition

[TABLE]

It is easy to verify from (2.15) that the scattered fields satisfy

[TABLE]

where $G_{\kappa_{j}}$ is the free space Green’s function for the two-dimensional Helmholtz equation at the wavenumber $\kappa_{j}$. Adding the incident field on both sides of (2.16a) and noting $\phi^{(2)}=\psi^{(2)}$ for (2.16b), we obtain

[TABLE]

Evaluating (2.17) at ${\bf r}_{i}$ and excluding the self-interaction, we get a nonlinear system of equations for $\phi_{k}^{(j)}$:

[TABLE]

which is the Foldy–Lax formulation for point scatterers with cubic nonlinearity.

3. Boundary integral formulation

In this section, we briefly introduce the boundary integral equation method for solving the scattering problem with extended scatterers. The detailed information can be found in [19].

Consider the scattering of a plane incident wave by a sound-soft extended scatterer, which is described by the domain $D$ with a boundary $\Gamma$ consisting of a finite number of disjoint, closed, bounded surfaces belonging to the class $C^{2}$. The exterior $\mathbb{R}^{2}\setminus\bar{D}$ is assumed to the connected, whereas $D$ itself is allowed to have more than one component. The unit normal vector $\nu$ on $\Gamma$ is assumed to be directed into the exterior of $D$.

The total field satisfies the Helmholtz equation:

[TABLE]

The sound-soft boundary implies that

[TABLE]

The total field $\phi$ consists of the incident field $\phi_{\rm inc}$ and the scattered field $\psi$:

[TABLE]

It follows from (2.2), (3.1), and (3.3) that the scattered field $\psi$ also satisfies the Helmholtz equation:

[TABLE]

The following Sommerfeld radiation condition is imposed to ensure the well-posedness of the scattering problem:

[TABLE]

It is shown in [19] by using the second Green theorem that

[TABLE]

and

[TABLE]

Adding the above two equations and using the sound-soft boundary condition (3.2), we get

[TABLE]

and

[TABLE]

To compute the scattered field $\psi$, it is required to compute $\partial_{\nu}\phi$ on $\Gamma$.

Taking the normal derivative of (3.6) on $\Gamma$ and applying the jump condition yield

[TABLE]

Multiplying (3.6) by ${\rm i}\eta$ and subtracting it from (3.8), we thus obtain a boundary integral equation for $\partial_{\nu}\phi$ on $\Gamma$:

[TABLE]

where the coupling parameter $\eta>0$ is introduced to ensure the unique solvability of (3.9).

4. Generalized Foldy–Lax formulation

This section presents the generalized Foldy–Lax formulation for the scattering problem of mixed scatterers, which consist of both the extended and point scatterers. We introduce the generalized Foldy–Lax formulation for the linear, quadratically nonlinear, and cubically nonlinear point scatterers, respectively.

4.1. Linear point scatterers

Viewing the external field acting on the point scatterers as point sources for the extended obstacle, we have the following equation for the total field:

[TABLE]

where $\phi_{k}$ is the external field acting on the $k$-th point scatterer and $\delta$ is the Dirac delta function. The obstacle is still assumed to be sound-soft. The total field vanishes on the boundary, i.e.,

[TABLE]

Subtracting the incident field (2.2) from the total field (4.1), we get the equation for the scattered field:

[TABLE]

As usual, the scattered field is required to satisfy the Sommerfeld radiation condition:

[TABLE]

We can follow the same steps as those in [19] to show that

[TABLE]

and

[TABLE]

Adding the above two equations and using the boundary condition (4.2), we have

[TABLE]

and

[TABLE]

To compute the scattered field $\psi$, it is required to compute $\partial_{\nu}\phi$ and $\phi_{k},k=1,\dots,m$.

Adding the incident field on both sides of (4.6) yields

[TABLE]

Evaluating (4.7) at ${\bf r}_{i}$ and excluding the self-interaction of the point scatterers, we get

[TABLE]

Taking the normal derivative of (4.5) on $\Gamma$ and applying the jump condition lead to

[TABLE]

Multiplying (4.5) by ${\rm i}\eta$ and subtracting it from (4.9), we obtain

[TABLE]

The coupled system (4.8) and (4.10) forms the generalized Foldy–Lax formulation for the scattering problem with the linear point scatterers and extended scatterers.

4.2. Quadratically nonlinear point scatterers

Consider the point scatterers with weak quadratic nonlinearity. Let $\kappa_{j}=\omega_{j}/c$ and $\phi^{(j)}$ be the field corresponding to the wavenumber $\kappa_{j}$. Viewing the external field acting on the nonlinear point scatterers as point sources for the extended obstacle, we have the following equations in the exterior domain $\mathbb{R}^{2}\setminus\bar{D}$:

[TABLE]

The sound-soft boundary condition implies that

[TABLE]

The fields satisfy the following relationship:

[TABLE]

where $\psi^{(j)}$ is the scattered field corresponding to the wavenumber $\kappa_{j}$ and satisfies the Sommerfeld radiation condition

[TABLE]

Similarly, we may show that the incident field satisfies

[TABLE]

the scattered fields satisfy

[TABLE]

and

[TABLE]

Adding the above equations and using the boundary condition yields

[TABLE]

and

[TABLE]

Adding the incident field to (4.14) and noting $\phi^{(2)}=\psi^{(2)}$, we obtain

[TABLE]

Evaluating (4.15) at ${\bf r}_{i}$ leads to

[TABLE]

Taking the normal derivative of (4.13) and using the jump conditions, we get

[TABLE]

Multiplying (4.13) by ${\rm i}\eta$ and subtract it from (4.17) give

[TABLE]

The coupled system (4.16) and (4.18) gives the generalized Foldy–Lax formulation for the scattering problem with quadratic nonlinear point scatterers and extended scatterers.

4.3. Cubically nonlinear point scatterers

Consider the point scatterers with weak cubic nonlinearity. Let $\kappa_{j}=\omega_{j}/c$ and $\phi^{(j)}$ be the field corresponding to the wavenumber $\kappa_{j}$. Viewing the external field acting on the nonlinear point scatterers as point sources for the extended obstacle, we have the following equations in the exterior domain $\mathbb{R}^{2}\setminus\bar{D}$:

[TABLE]

The sound-soft boundary condition implies that

[TABLE]

The fields satisfy the following relationship:

[TABLE]

where $\psi^{(j)}$ is the scattered field corresponding to the wavenumber $\kappa_{j}$ and satisfies the Sommerfeld radiation condition

[TABLE]

Similarly, we may show that the incident field satisfies

[TABLE]

the scattered fields satisfy

[TABLE]

and

[TABLE]

Adding the above equations and using the boundary condition yield

[TABLE]

and

[TABLE]

Adding the incident field to (4.22) and noting $\phi^{(3)}=\psi^{(3)}$, we obtain

[TABLE]

Evaluating (4.23) at ${\bf r}_{i}$ and excluding the self-interaction lead to

[TABLE]

Taking the normal derivative of (4.21) and using the jump conditions, we get

[TABLE]

Multiplying (4.13) by ${\rm i}\eta$ and subtract it from (4.17) give

[TABLE]

The coupled system (4.24) and (4.26) gives the generalized Foldy–Lax formulation for the scattering problem with cubic nonlinear point scatterers and extended scatterers.

5. Direct imaging method

In this section, we introduce a fast direct imaging method to reconstruct the shape of the extended scatterers.

5.1. Far-field pattern

The far-field pattern is an important quantity which encodes the information about the scatterers such as location and shape. Given an incident field with incident direction ${\bf d}$, the scattered field has the asymptotic behavior

[TABLE]

uniformly in all directions $\hat{\bf r}={\bf r}/|{\bf r}|$, where the function $\psi_{\infty}$ is called the far-field pattern of the scattered field $\psi$, and $\hat{\bf r}\in\mathbb{S}^{1}$ is the observation direction.

Recall the asymptotic behavior for the Hankel function for large arguments

[TABLE]

and the following identity

[TABLE]

Using (5.1) and the scattered field representations (2.6), (2.11), (2.16), (3.7), (4.6), (4.14), (4.22), we obtain the following far-field patterns of the scattered field for the scattering problem with point scatterers, extended scatterers, and mixed scatterers, respectively.

1. (i)

Foldy–Lax formulation for point scatterers:

1. (a)

Linear point scatterers

[TABLE] 2. (b)

Quadratically nonlinear point scatterers

[TABLE] 3. (c)

Cubically nonlinear point scatterers

[TABLE] 2. (ii)

Boundary integral formulation for extended scatterers

[TABLE] 3. (iii)

Generalized Foldy–Lax formulation for mixed scatterers

1. (a)

Linear point scatterers

[TABLE] 2. (b)

Quadratically nonlinear point scatterers

[TABLE] 3. (c)

Cubically nonlinear point scatterers

[TABLE]

where

[TABLE]

When the observation directions and the number of point scatterers are large, it is very slow to directly evaluate the far-field patterns (5.2)–(5.8). In practice, the evaluation of the far-field patterns are accelerated by the fast multipole method (FMM)[26].

5.2. Imaging function

Consider an array of transmitters that can send out plane incident waves and record the far-field pattern of the scattered waves. Assume that we have a set of incident plane waves with incident directions ${\bf d}_{1},\dots,{\bf d}_{N}$ and the far-field patterns are recorded at observation direction $\hat{\bf r}_{1},\dots,\hat{\bf r}_{M}$, where ${\bf r}_{i}=(\cos\alpha_{i},\sin\alpha_{i})$ and ${\bf d}_{j}=(\cos\beta_{j},\sin\beta_{j}),i=1,\dots,M,j=1,\dots,N$. Here $\alpha_{i}$ is the observation angle and $\beta_{j}$ is the incident angle. These measurement of the far-field patterns form an $M\times N$ response matrix

[TABLE]

where the far-field pattern $\psi_{\infty}$ represents any one of the far-field patterns in (5.2)–(5.8).

Consider two unit vectors:

[TABLE]

and

[TABLE]

which is the illumination vector with respect to the receivers and the transmitters, respectively. Define an imaging function

[TABLE]

The direct imaging method is to evaluate the imaging function (5.10) at any given sampling point ${\bf r}\in\mathbb{R}^{2}$.

Since the imaging function $I({\bf r})$ needs to be evaluated at every sampling point and each evaluation requires the matrix-vector multiplication, the computation is intensive. However, the imaging function can be written as

[TABLE]

It is clear to note from (5.2) that the imaging function $I_{\kappa}({\bf r})$ is the two-dimensional Fourier transform of the $(i,j)$-th entry of the response matrix $P_{\kappa}^{(i,j)}$, which takes the response matrix from the frequency space $(\kappa(\cos{\alpha_{i}}+\cos{\beta_{j}}),\kappa(\sin{\alpha_{i}}+\sin{\beta_{j}}))$ to the physical space $(x,y)$. While the frequency points $(\kappa(\cos{\alpha_{i}}+\cos{\beta_{j}}),\kappa(\sin{\alpha_{i}}+\sin{\beta_{j}}))$ are not uniform, one needs to apply the two-dimensional non-uniform fast Fourier transform (NUFFT) to accelerate the evaluation.

5.3. NUFFT

Here we give a short introduction to NUFFT in one dimension. Details can be found in [25]. Application of NUFFT on the acceleration of evaluating forward multi-particle scattering in a layered medium can be found in[36].

Consider the following expression:

[TABLE]

where $c_{j}\in\mathbb{C},x_{i}\in[-m/2,m/2-1],\xi_{j}\in[-\pi,\pi]$. Given $c_{j}$ and $\xi_{j}$, the goal is to evaluate $f(x_{i})$ efficiently. When $m=n$ and $x_{i}$, $\xi_{j}$ are uniformly distributed in their intervals, the sum can be accelerated by standard FFT.

Consider the situation where $m\neq n$ and $\xi_{j}$ is not uniformly given in the interval $[-\pi,\pi]$. Assume that $x_{i}$ is still uniformly distributed. We first rewrite equation (5.12) as

[TABLE]

where $\delta$ is the Dirac delta function and

[TABLE]

In other words, $f(x)$ can be taken as the Fourier transform of $F(\xi)$. However, $F(\xi)$ consists of delta functions, which is numerically difficult to evaluate. The function $F(\xi)$ can be modified by convolving it with Gaussian function $g(\xi)=e^{-\xi^{2}/{4\tau}}$, where $\tau>0$ is a small number. Let

[TABLE]

We over sample $\tilde{F}(\xi)$ uniformly on the interval $[-\pi,\pi]$ and take the Fourier transform of $\tilde{F}(\xi)$ with the standard FFT. Notice that to apply FFT, both $F(\xi)$ and $g(\xi)$ are periodized first on $[-\pi,\pi]$ before evaluating $\tilde{F}(\xi)$[25]. Denote by $\hat{\tilde{F}}(x)$ and $\hat{g}(x)$ the Fourier transform $\tilde{F}(\xi)$ and $g(\xi)$, respectively. It follows from the convolution property of Fourier transform that $f(x)$ can be approximated by

[TABLE]

To conclude, the NUFFT consists of three steps: (1) convolution with Gaussian function, (2) standard FFT for the over sampled function, (3) deconvolution with Gaussian function. The choices of parameter $\tau$ and oversampling factor in step 2 highly affect the numerical performance of the NUFFT. For the one dimensional case, $\tau$ is chosen to be $12/m^{2}$ and the number of over sampling points is $2m$, which guarantees 12 digits precision. Using these parameters, we can easily find out that the complexity of the one-dimensional NUFFT is $O(K\log K)$, where $K=\max\{m,n\}$. The scheme introduced here can be easily extended to the two- and higher-dimensional cases. The two-dimensional NUFFT is used in this work to evaluate (5.2).

5.4. Imaging with nonlinear point scatterers

The method introduced in the beginning of this section works effectively to image the extended scatterers surrounded by linear point scatterers. The method can also be applied directly to the scattered field generated by the first order frequency wave in the nonlinear case. The purpose of adding nonlinear point scatterers is to introduce higher order field, which allows to capture more details of the extended scatterers.

However, if we fix the location of the nonlinear point scatterers, the imaging method will fail for testing the higher order frequency waves. The reason is obvious: in the direct imaging method, the test function ${\bf v}_{\kappa}$ depends on the incident angle ${\bf d}_{j},j=1,\dots,N$. For higher order frequency waves, the incident wave is essentially generated by the nonlinear interaction of point scatterers. The direction of such generated incident wave does not line up with the incident direction in the test function.

To remedy the scheme, we first put the point scatterers sufficiently far away from the extended scatterers and line up the point scatterers towards the incident direction, as is shown in Figure 2. In addition, we move the point scatterers along with the change of incident direction. For the higher order frequency waves, the entries of the response matrix are taken as the difference of the far-field patterns from the generalized Foldy–Lax formulation and the Foldy–Lax formulation, i.e.,

[TABLE]

for the quadratically nonlinear point scatterers and

[TABLE]

for the cubically nonlinear point scatterers. The purpose of taking the difference is to avoid the possibility that the far-field pattern from the point scatterers may dominate the far-field pattern from the extended scatterers.

6. Numerical experiments

In this section, we discuss the implementation of the direct scattering problem and present some numerical experiments for the inverse scattering problem. In all of the following examples, the extended scatterer is a five-leaf shaped obstacle and can be parameterized, up to a shift and rotation, by

[TABLE]

where $t\in[0,2\pi]$ is the parameter. For convenience, we summarize some of the parameters used in the numerical experiments in Table 1. The resulting system of equations are obtained after discretizing the boundary of the extended scatterers. The total number of unknowns is $N_{\rm point}+N_{\rm boundary}$ for linear point scatterers and $2(N_{\rm point}+N_{\rm boundary})$ for nonlinear point scatterers.

6.1. Direct scattering problem solver

Let $\boldsymbol{\phi}=(\phi_{1},\dots,\phi_{m})^{\top},\boldsymbol{\phi}_{\rm inc}=(\phi_{\rm inc}({\bf r}_{1}),\dots,\phi_{\rm inc}({\bf r}_{m}))^{\top}$. Define an $m\times m$ matrix

[TABLE]

and three linear operators

[TABLE]

The generalized Foldy–Lax formulation of linear point scatterers (4.8) and (4.10) can be written as the operator form:

[TABLE]

where $\varphi=\partial_{\nu}\phi({\bf r})$ and $\varphi_{\rm inc}=(\partial_{\nu}-{\rm i}\eta)\phi_{\rm inc}({\bf r})$. In the discretization, the surface of the extended obstacle $\Gamma$ is discretized by a set of uniform points in the parameter space $t$; the singular integral is evaluated by the Alpert quadrature [4]; the boundary integral equations are solved by the Nyström method.

There are three approaches to solve the linear system (6.2):

1. (1)

Direct solver. Apply the LU factorization from the Lapack library and parallelize it by OpenMP on a multicore workstation; 2. (2)

Iterative solver. Apply GMRES to the whole system and accelerate the matrix vector product by the fast multipole method (FMM)[26]. 3. (3)

Hybrid method. Assume that the number of points to discretize the extended obstacle is relatively small compared to the number of point scatterers. First is to invert $\mathcal{K}_{\kappa}$ directly and then solve the Schur complement of (6.2) by an iterative method, i.e., solve iteratively with the FMM acceleration of the linear system

[TABLE]

To investigate the three different methods, we solve the scattering problem for the extended obstacle in (6.1), which is surrounded by a group of linear point scatterers. Table 2 shows the numerical performance for the three different methods. Obviously, the direct solver is the best option for the linear problem. It solves the system very rapidly with parallelization, and it is independent of the location of the point scatterers and the wavenumber $\kappa$. Both of the iterative methods fail if the point scatterers are randomly and densely distributed in a specific area. In addition, to solve the inverse problem, the direct problem (6.2) has to be solved many times with different right hand sides in order to construct the response matrix, which corresponds to the far-field patterns for different incident directions. The advantage of the direct solver is clear: we only need to invert (6.2) once and apply matrix vector product for different right hand sides, which can also be done by parallelization.

Next we discuss the solver for the generalized Foldy–Lax of nonlinear point scatterers. We only describe the steps for the quadratically nonlinear point scatterers since they are similar to the cubically nonlinear point scatterers.

Let $\boldsymbol{\phi}^{(j)}=(\phi_{1}^{(j)},\dots,\phi_{m}^{(j)})^{\top}$ and $\varphi^{(j)}=\partial_{\nu}\phi^{(j)}({\bf r})$ be the external field acting on the point scatterers and the normal derivative of the total field on the boundary of the extended obstacle at the wavenumber $\kappa_{j}$, respectively. The generalized Foldy–Lax for nonlinear point scatterers can be written as

[TABLE]

where $\mathcal{D}_{ij}$ represents the nonlinear interaction between $\boldsymbol{\phi}^{(1)}$ and $\boldsymbol{\phi}^{(2)}$ at the point scatterers, $\mathcal{H}_{ij}$ is the nonlinear interaction from the point scatterers to the extended obstacle, $\mathcal{M}_{\kappa_{j}}$ denotes the linear interaction from the extended obstacle to the point scatterers at the wavenumber $\kappa_{j}$, and $\mathcal{K}_{\kappa_{j}}$ is the linear interaction for the extended obstacle.

Due to the large number of unknowns and nonlinearity, neither the direct solver nor the iterative method is applicable to the nonlinear system (6.3). Assuming that the number of point scatterers is relatively small, we propose an efficient nonlinear solver, which can be be applied to the Schur complement of (6.3). Specifically, we invert $\mathcal{K}_{\kappa_{1}}$ and $\mathcal{K}_{\kappa_{2}}$ directly, and then solve the following nonlinear system:

[TABLE]

Finally, we apply a trust region Newton type method to solve the resulting nonlinear system (6.1).

6.2. Linear point scatterers

We show the imaging results of extended scatterers which are surrounded by a group of linear point scatterers. These point scatterers are randomly and uniformly distributed in the annulus $\{{\bf r}\in\mathbb{R}^{2}:10<|{\bf r}|<11\}$. The scattering coefficients are $\sigma_{k}=0.5$ for all the point scatterers. The numerical performance of the following two examples is shown in Table 3. It is clear to note the the proposed method is not only efficient for the direct problem simulation but also for the inverse problem imaging.

6.2.1. Example 1

This example is to image two extended scatterers which are surrounded by 1000 linear point scatterers at the wavenumber $\kappa=10$. The imaging result is shown in Figure 3(a). The extended scatterers are well reconstructed except for the parts where these two scatterers are close.

6.2.2. Example 2

To show the influence of the wavenumber on the imaging resolution, we take the same extended and point scatterers as those in Example 1, but we use the wavenumber $\kappa=50$. The imaging result is shown in Figure 3(b), which has a better resolution than Figure 3(a) does. These two scatterers are well reconstructed even for the parts where they are close. As is expected, higher wavenumber can capture finer structures.

6.3. Quadratically nonlinear point scatterers

We consider the imaging of extended scatterers with two quadratically nonlinear point scatterers. The nonlinear scattering coefficients are $\sigma_{k,1}^{(1)}=\sigma_{k,2}^{(1)}=0.5$ and $\sigma_{k,1}^{(2)}=\sigma_{k,2}^{(2)}=0.4$. The numerical performance is shown in Table 4 for the following two examples. As can be seen, the proposed method is also efficient for the nonlinear direct problem simulation.

6.3.1. Example 3

This example is to image one extended scatterers with two quadratically nonlinear point scatterers. The wavenumber of the incident wave is $\kappa=2$. First we consider the case when the two point scatterers are fixed at the location $(-13,0)$ and $(-14,0)$, respectively. The imaging result is shown in Figure 4. The linear wave can reconstruct the extended scatterer but with poor resolution. The wave from the second harmonic generation cannot reconstruct the extended scatterer. As is described in section 4.4, we change the location of the two point scatterers on the circle with radius 13 and 14 aligned with the incident direction. For instance, if the angle of the incidence is $\theta=\frac{\pi}{3}$, the two point scatterers are located at $(13\cos\theta,13\sin\theta)$ and $(14\cos\theta,14\sin\theta)$. The imaging result is shown in Figure 5. The linear wave yields almost the same imaging result as the fixed point scatterers does. But the wave from the second harmonic generation gives a much better imaging result due to the doubled frequency.

6.3.2. Example 4

This example is to image two extended scatterers with two point scatterers located on the circles of radii $|{\bf r}|=13$ and $|{\bf r}|=14$. The wavenumber of the incident wave is $\kappa=5$. The imaging result is shown in Figure 6. The linear wave can still produces a reasonable imaging result. But the nonlinear wave fails to identify the two extended scatterers. Next we keep everything else the same except moving away the two point scatterers to the circles of radii $|{\bf r}|=130$ and $|{\bf r}|=131$. The imaging result is shown in Figure 7. It can be seen that the nonlinear wave can clearly identify the two extended scatterers. The reason for this is clear: when the point scatterers are far away, their generated wave interacted with the extended obstacles is more like a plane wave incidence, which can be better resolved by the illumination vectors.

6.4. Cubically nonlinear point scatterers

Finally, we investigate the performance of using cubically nonlinear point scatterers. We put two cubically nonlinear point scatterers around the extended scatterers. The scattering coefficients are $\sigma_{k,1}^{(1)}=\sigma_{k,2}^{(1)}=0.5$ and $\sigma_{k,1}^{(3)}=\sigma_{k,2}^{(3)}=\sigma_{k,3}^{(3)}=0.4$. The numerical performance is summarized in Table 5.

6.4.1. Example 5

This example is to image one extended scatterer with two point scatterers located at the circles with radii $|{\bf r}|=13$ and $|{\bf r}|=14$, respectively. The wavenumber of the incidence is $\kappa=2$. The imaging result is shown in Figure 8. Comparing Figure (5)(b) with Figure (8)(b), we observe that the cubically nonlinear point scatterers can produce a better resolution than the quadratically nonlinear point scatterers does, which confirms once again that higher frequency produces finer resolution.

6.4.2. Example 6

This example is to image two extended scatterers with two point scatterers located at the circles with radii $|{\bf r}|=13$ and $|{\bf r}|=14$, respectively. The wavenumber of the incidence is $\kappa=5$. The imaging result is shown in Figure 9. Comparing Figure (7)(b) with Figure (9)(b), we observe the same pattern that higher frequency wave generates better resolved imaging result.

7. Conclusion

We have presented a generalized Foldy–Lax formulation for the scattering by the heterogeneous media consisting of linear or nonlinear point scatterers and extended obstacles. A new imaging function is proposed and a fast direct imaging method is developed for the inverse obstacle scattering problem. Using the nonlinear point scatterers to excite high harmonic generation, enhanced imaging resolution is achieved to reconstruct the extended obstacle. Our method shares the attractive features of all the other direct imaging methods. In addition, the evaluation is accelerated by using the FFT due to the special construction of the imaging function. Numerical results show that the method is effective to solve the inverse obstacle scattering problem. The proposed method can be directly extended to solve the inverse obstacle scattering problem in three dimensions. As an FFT-based method, a greater potential can show up in higher dimensions. We are currently working on the three-dimensional problem and the results will be reported elsewhere.

Appendix A Nonlinear wave equations

In the scalar theory of electromagnetic fields, the scalar electric field $u(\textbf{r},t)$ obeys the wave equation:

[TABLE]

where $c>0$ is the light speed and $P({\bf r},t)$ is the polarization density.

We adopt the following Fourier and inverse Fourier transformation convention:

[TABLE]

Note that if $f(\textbf{r},t)$ is a real-valued function, then $f(\textbf{r},-\omega)=\bar{f}(\textbf{r},\omega)$, where the bar is the complex conjugate. Taking the Fourier transform of (A.1), we obtain the Helmholtz equation:

[TABLE]

where $\kappa(\omega)=\omega/c$ is the wavenumber.

The polarization may be expanded in powers of the electric field. In principle the expansion involves infinitely many terms, but only the first few terms are of practical importance if the nonlinearity is weak. In this paper, we consider linear, quadratically nonlinear, and cubically nonlinear media.

(i) A medium is linear if

[TABLE]

where the coefficient $\chi^{(1)}(\textbf{r},\omega)$ is the first-order susceptibility. Combining (A.2) and (A.3) gives

[TABLE]

(ii) A medium is quadratically nonlinear if

[TABLE]

where $\chi^{(2)}(\textbf{r},\omega_{1},\omega_{2})$ are the second-order susceptibilities. The summation indicates that the electric fields at the frequencies $\omega_{1}$ and $\omega_{2}$ contribute to the polarization at the frequency $\omega$ if $\omega_{1}+\omega_{2}=\omega$. Second-order nonlinear effects include second harmonic generation, which is excited by a monochromatic incident field of frequency $\omega$ in a quadratically nonlinear medium. We assume that the nonlinear susceptibilities are sufficiently weak, i.e.,

[TABLE]

We also assume that the second-order susceptibilities have full permutation symmetry, i.e.,

[TABLE]

Let $\omega_{1}=\omega$ and $\omega_{2}=2\omega$. Define $\kappa_{j}=\omega_{j}/c$. Denote by $u^{(j)}$ the field corresponding to the wavenumber $\kappa_{j}$. It follows from (A.2) and (A.5) that $u^{(j)}$ satisfies

[TABLE]

(iii) A medium is cubically nonlinear if

[TABLE]

where $\chi^{(3)}(\textbf{r},\omega_{1},\omega_{2},\omega_{3})$ are the third order susceptibilities. Materials with inversion symmetry have zero second order susceptibilities and thus fall into this category. Third-order nonlinear effects include third-harmonic generation. We assume that the nonlinear susceptibilities are sufficiently weak, i.e.,

[TABLE]

In addition, we assume that third-order susceptibilities have full permutation symmetry, i.e.,

[TABLE]

where $\{{\rm p}(1),{\rm p}(2),{\rm p}(3)\}$ is a permutation of $\{1,2,3\}$. Suppose that a source of frequency $\omega$ is incident upon a cubic nonlinear medium. Let $\omega_{1}=\omega$ and $\omega_{3}=3\omega$. Define $\kappa_{j}=\omega_{j}/c$. Denote by $u^{(j)}$ the field corresponding to the wavenumber $\kappa_{j}$. It follows from (A.2) and (A.7) that $u^{(j)}$ satisfies

[TABLE]

Appendix B Point scatterer models

Scatterers that are small compared to the wavelength can be effectively treated as point scatterers [20]. The susceptibilities of a small scatterer located at ${\bf r}_{0}$ can be replaced by delta functions in the following forms:

[TABLE]

where $\eta^{(j)}$ are the effective susceptibilities of the point scatterer. The representations (B.1) are valid outside of the interaction region, which is of size comparable to the physical size of the scatterer.

(i) For a collection of linear point scatterers located at ${\bf r}_{k}\in\mathbb{R}^{2},k=1,\dots,m$, the wave equation (A.4) takes the form:

[TABLE]

Identifying $u(\textbf{r},\omega)$ with $\phi(\bf r)$ and $4\pi\kappa^{2}(\omega)\eta^{(1)}({\bf r}_{k},\omega)$ with $\sigma_{k}$, we obtain (2.3).

(ii) For a collection of quadratically nonlinear point scatterers located at ${\bf r}_{k}\in\mathbb{R}^{2},k=1,\dots,m$, the wave equation (A.6) takes the form:

[TABLE]

Identifying $u^{(j)}$ with $\phi^{(j)}$, $4\pi\kappa_{1}^{2}\eta^{(1)}({\bf r}_{k},\omega_{1})$ with $\sigma^{(1)}_{k,1}$, and $8\pi\kappa_{1}^{2}\eta^{(2)}({\bf r}_{k},\omega_{2},-\omega_{1})$ with $\sigma^{(2)}_{k,1}$, we obtain (2.9a). Identifying $4\pi\kappa_{2}^{2}\eta^{(1)}({\bf r}_{k},\omega_{2})$ with $\sigma^{(1)}_{k,2}$, and $4\pi\kappa_{2}^{2}\eta^{(2)}({\bf r}_{k},\omega_{1},\omega_{1})$ with $\sigma^{(2)}_{2,k}$, we obtain (2.9b). Note that $\sigma^{(1)}_{k,1}$ and $\sigma^{(1)}_{k,2}$ generically take different values, so do $\sigma^{(2)}_{k,1}$ and $\sigma^{(2)}_{k,2}$.

(iii) For a collection of cubically nonlinear point scatterers located at ${\bf r}_{k}\in\mathbb{R}^{2},k=1,\dots,m$, the wave equation (A.8) takes the form:

[TABLE]

Identifying $u^{(j)}$ with $\phi^{(j)}$, $4\pi\kappa_{1}^{2}\eta^{(1)}({\bf r}_{k},\omega_{1})$ with $\sigma^{(1)}_{k,1}$, $12\pi\kappa_{1}^{2}\eta^{(3)}({\bf r}_{k},\omega_{1},\omega_{1},-\omega_{1})$ with $\sigma^{(3)}_{k,1}$, and $12\pi\kappa_{1}^{2}\eta^{(3)}({\bf r}_{k},\omega_{3},-\omega_{1},-\omega_{1})$ with the second $\sigma^{(3)}_{k,2}$, we obtain (2.14a). Identifying and $4\pi\kappa_{3}^{2}\eta^{(1)}({\bf r}_{k},\omega_{3})$ with $\sigma^{(1)}_{k,2}$ and $4\pi\kappa_{3}^{2}\eta^{(3)}({\bf r}_{k},\omega_{1},\omega_{1},\omega_{1})$ with $\sigma^{(3)}_{k,3}$ , we obtain (2.14b). Note that $\sigma^{(1)}_{k,1}$ and $\sigma^{(1)}_{k,2}$ generically take different values, and so do $\sigma^{(1)}_{k,1}$, $\sigma^{(1)}_{k,2}$ and $\sigma^{(1)}_{k,3}$.