## Abstract

The formal representation of the quasi-periodic Helmholtz Green's function obtained by the method of images is only conditionally convergent and, thus, requires an appropriate summation convention for its evaluation. Instead of using this formal sum, we derive a candidate Green's function as a sum of two rapidly convergent series, one to be applied in the spatial domain and the other in the Fourier domain (as in Ewald's method). We prove that this representation of Green's function satisfies the Helmholtz equation with the quasi-periodic condition and, furthermore, leads to a fast algorithm for its application as an operator.

We approximate the spatial series by a short sum of separable functions given by Gaussians in each variable. For the series in the Fourier domain, we exploit the exponential decay of its terms to truncate it. We use fast and accurate algorithms for convolving functions with this approximation of the quasi-periodic Green's function. The resulting method yields a fast solver for the Helmholtz equation with the quasi-periodic boundary condition. The algorithm is adaptive in the spatial domain and its performance does not significantly deteriorate when Green's function is applied to discontinuous functions or potentials with singularities. We also construct Helmholtz Green's functions with Dirichlet, Neumann or mixed boundary conditions on simple domains and use a modification of the fast algorithm for the quasi-periodic Green's function to apply them.

The complexity, in dimension *d*≥2, of these algorithms is (*κ*^{d} log *κ*+*C*(log *ϵ*^{−1})^{d}), where *ϵ* is the desired accuracy, *κ* is proportional to the number of wavelengths contained in the computational domain and *C* is a constant. We illustrate our approach with examples.

## 1. Introduction

In this paper, we construct fast and accurate algorithms for applying Helmholtz Green's functions incorporating a variety of boundary conditions on simple domains. We come to this problem from the perspective of developing fast algorithms for applying operators of mathematical physics with finite but arbitrary accuracy (in operator norm) as we discuss later. Instead of emphasizing computing of values of Green's functions, as is typical in the literature mentioned below, we focus on the problem of applying such Green's functions as operators in a fast and accurate manner. We note that the accurate computation of the values of a Green's function does not by itself resolve the issue of its efficient application and use as an operator. Towards this end, we develop approximations of Green's functions that resolve the problem of algorithmic efficiency in applying them to discontinuous functions or potentials with singularities.

The key element of our approach is a fast algorithm for computing convolutions with the quasi-periodic Helmholtz Green's function,(1.1)for functions , where *D* is the primitive cell of a Bravais lattice *Λ*. We consider the case where the lattice is defined by *d* linearly independent vectors in dimension *d*≥2. Green's function *G*_{q} satisfies(1.2)(1.3)where *κ*>0, ** l**∈

*Λ*,

**∈**

*x**D*and

**∈**

*k*^{d}. The quasi-periodicity vector

**is sometimes referred to as Bloch or crystal momentum vector. Consequently, the function**

*k**u*in (1.1) satisfies(1.4)(1.5)It is well known that the method of images (typically used in dimensions

*d*=2 or 3) gives rise to a natural, yet formal representation,(1.6)where |.| denotes the length of a vector and is the zeroth-order Hankel function of the first kind. These sums are only conditionally convergent; thus, they require a summation convention to yield a classical solution of (1.4) and (1.5).

The formal sum (1.6) appears in many areas of physics and engineering. For example, one of the first applications (Ewald (1913), see Cruickshank *et al*. (1992) for an English translation) describes X-ray diffraction by crystals using (1.6) as the key mathematical object. In fact, by transforming (1.6) to the Fourier domain, it is easy to obtain directions and X-ray frequencies by using the so-called Ewald's sphere. The sum (1.6) occurs in wave propagation in periodic structures (e.g. Brillouin 1953), in the study of photonic crystals (e.g. Soukoulis 1992; Joannopoulos *et al*. 1995), in band structure computations in solid-state physics (e.g. Kohn & Rostoker 1954; Ham & Segall 1961), as well as in the analysis of ergodic systems (Berry 1981).

Besides Ewald's (1921) method, there are other approaches for interpreting and evaluating (1.6) (see Glasser & Zucker (1980) and Linton (1998) for a survey and references therein). One approach is to use addition theorems for Bessel functions, which separates the input variables ** x** from the lattice vectors

**. In such cases, Green's function is written as an expansion with respect to Bessel functions with coefficients given by a lattice sum. We note that Yasumoto & Yoshitomi (1999), McPhedran**

*l**et al*. (2000) and Dienstfrey

*et al*. (2001) employ this approach and develop different algorithms for evaluating the lattice sum. Although these algorithms yield a method for the fast evaluation of Green's function, they do not provide an algorithm to perform fast convolutions.

In our approach, as in Ewald's (1921) summation, we split Green's function into two absolutely convergent series, one in the spatial domain and the other in the Fourier domain. We then verify directly that our representation of Green's function satisfies (1.2) and (1.3) in the usual sense and, thus, (1.1) yields the classical solution of (1.4) and (1.5). We also provide an alternative derivation of Ewald's method using a limiting procedure based on analytic continuation, which more naturally connects with our approach and goals.

Using the periodic Green's function (** k**=

**0**in (1.3)) and the method of images, we construct other Helmholtz Green's functions that incorporate either Dirichlet, Neumann or mixed boundary conditions on simple domains. The resulting integral operators are no longer convolutions, but the algorithm for applying these Green's functions is similar to that for the quasi-periodic Green's function. The application of Green's functions satisfying boundary conditions is also split between the spatial and the Fourier domains resulting in an algorithm with the same computational complexity as that for the quasi-periodic Green's function.

In the spatial domain, we approximate operators using separated representations given by a sum of Gaussians and note that this type of approximation (e.g. Beylkin & Mohlenkamp 2002, 2005; Beylkin & Monzón 2005) has been successfully used by Harrison *et al*. (2004) and Beylkin *et al*. (2007, 2008) to construct fast and accurate algorithms for applying non-oscillatory kernels. In this paper, we extend the results in Beylkin *et al*. (submitted) to the quasi-periodic Helmholtz Green's functions as well as Green's functions with boundary conditions on simple domains. Thus, to apply Green's function, in the spatial domain we compute convolutions with Gaussians via an adaptive multiresolution algorithm (e.g. Harrison *et al*. 2004; Beylkin *et al*. 2007, 2008), or the fast Gauss transform in Greengard & Strain (1991), Strain (1991) and Greengard & Sun (1998).

In the Fourier domain, we use the unequally spaced fast Fourier transform (USFFT) (e.g. Dutt & Rokhlin 1993; Beylkin 1995; Lee & Greengard 2005). The resulting algorithm computes (1.1) with controlled accuracy *ϵ*, computational cost proportional to (log *ϵ*^{−1})^{d}, and maintains its performance when applied to functions with discontinuities or singularities (e.g. Coulomb or Lennard-Jones potentials). The same approach yields algorithms for Green's functions in simple domains with either Dirichlet, Neumann or mixed boundary conditions. These algorithms yield fast and adaptive solvers for the Helmholtz equation with the aforementioned boundary conditions and have computational complexity , where *C* is a constant.

We proceed by providing some preliminaries in §2 and, in §3, derive a representation of the quasi-periodic Green's function as a sum of two convergent series solving (1.2) and (1.3). Using this representation, in §4 we develop a fast algorithm for computing convolutions with the quasi-periodic Green's function and illustrate our approach with examples. Then, in §5, we construct approximations to Green's functions satisfying boundary conditions on simple domains. Finally, we discuss implications of our approach for developing a unified methodology for applying oscillatory and non-oscillatory Green's functions.

## 2. Preliminaries

### (a) Bravais lattice

A Bravais lattice in dimension *d*=2 is defined aswhere the lattice vectors are linearly independent. The reciprocal lattice (in the Fourier domain) is then given bywhere are the reciprocal lattice vectors defined to satisfySimilarly, in dimension *d*=3 we haveandwhere and . An obvious generalization yields Bravais lattices in any dimension *d*.

In dimension *d*, we consider the primitive cell to be the *d*-dimensional parallelepiped associated with the vectors . We denote the primitive cell as *D* and its volume as *V*. We refer to, for example, Kittel (1986) for a detailed description of Bravais lattices.

### (b) Fourier transform and Poisson summation formula

We use the Fourier transform in dimension *d*and its inverseFor our purposes, it is sufficient to consider the Schwartz class of functions (^{d}) containing infinitely differentiable functions with derivatives decaying faster than any inverse polynomial (e.g. Grafakos 2004, §2*b*). We have

(*Poisson summation formula*) *Let f*∈(^{d})*,* *Λ be a Bravais lattice, Λ*^{*} *the reciprocal lattice and V the volume of the primitive cell. Then**for* .

This result follows from observing that the set of functions forms an orthonormal basis for (*D*) since the linear change of variables , where , reduces the problem to that of the standard Fourier series .

### (c) Free-space Green's function

The outgoing free-space Helmholtz Green's function in dimension *d* (where is the Hankel function of the first kind),satisfies the Helmholtz equation(2.1)and the Sommerfeld radiation conditionIn particular, we haveOn taking the Fourier transform of (2.1), we obtainThe inverse Fourier transform of is a singular integral and we use the limiting procedure in Beylkin *et al*. (submitted) to define Green's function as(2.2)

## 3. Quasi-periodic Green's function via absolutely convergent series

The quasi-periodic Green's function formally described by (1.6) requires a summation convention for its evaluation. Instead, we construct the quasi-periodic Green's function via a sum of two convergent series that yield an explicit (classical) solution of (1.2) and (1.3) and, in §4, describe a fast algorithm for its application.

As a motivation, let us considerand note that this sum formally satisfies (1.2) and (1.3) provided , where . Following the approach in Beylkin *et al*. (submitted), we choose *η*>0 and split the above expression asWe show below that the splitting parameter *η* is the same as in Ewald's method and we discuss practical considerations for its selection in §4.

Observing that the first term is an absolutely convergent series, we define(3.1)for *κ*>0, , where . Note that the sum in (3.1) is independent of the dimension *d*.

For the second term, we use the Poisson summation and (still formally) obtain(3.2)where(3.3)We note that by replacing e^{−2s} by its maximum on [log (2*η*), ∞) and changing variables *t*=e^{2s}/4−*η*^{2}, we obtain the estimate with some constant *C* and any ** l**≠

**0**.

This estimate shows that (3.2) is an absolutely convergent series which we use as a definition for *G*_{spatial}. We define(3.4)as a sum of two absolutely convergent series with *G*_{spatial} in (3.2) and *G*_{Fourier} in (3.1).

Next, we show that *G*_{q} in (3.4) satisfies (1.2) and (1.3) in the usual sense and, therefore, we no longer depend on an interpretation of (1.6) as an operator. As a result, we may consider convolving *G*_{q} with functions from various classes, e.g. *L*^{p}(*D*), and the convolution (1.1) gives us a classical solution of (1.4) and (1.5). We prove that

*The function G*_{q} *in* (*3.4*) *is the quasi-periodic Green's function satisfying* (*1.2*) *and* (*1.3*) *for κ*>0*,* *,* *where* *and* . *This result holds in any dimension d*≥2.

The quasi-periodic condition for *G*_{Fourier} in (3.1) follows fromsince for any ** l**∈

*Λ*and . For

*G*

_{spatial}in (3.2), a shift in summation yields a factor e

^{−ik.l}and, thus, it also satisfies the quasi-periodic condition.

We apply Δ+*κ*^{2} to (3.1) and (3.2). Using proposition 2.1, we have(3.5)For *G*_{spatial} in (3.2), we change the variables in the integral, *s*=log(2*t*), and obtain(3.6)In the previous sum, we separate the ** l**=0 term and note that for we have(3.7)This identity follows by observing thatand using integration by parts. Thus, for

**∈**

*x**D*and

**≠0, we have a term-by-term cancellation between (3.6) and (3.5), which yieldswhere(3.8)The function**

*l**F*corresponds to the difference of

**=0 terms in (3.6) and (3.5).**

*l*Note that due to (3.7), *F*(** x**) vanishes for

**≠0. In order to show that**

*x**G*

_{q}satisfies (1.2), we show that the Fourier transform of (3.8) is . We haveandThus, we obtain ▪

In the special case where for some and fixed , where is the reciprocal lattice, each of the functions satisfies the quasi-periodic condition (1.5) and is a eigenfunction of Δ with eigenvalue . Thus, such functions are in the null space of Δ+*κ*^{2} and we require the function *f* in (1.1) to satisfy(3.9)for all such vectors ** d**. Therefore, the solution to (1.4) and (1.5) is given by(3.10)where

*c*

_{d}are arbitrary constants and

We note that to derive *G*_{q}, it is sufficient to consider the real part of the free-space Green's function *G*_{free} since the imaginary part of *G*_{free} leads to a vanishing sum. Indeed, using the Poisson summation for generalized functions, for *κ*>0, , and ** x**∈

*D*, we havewhere

*Λ*is a Bravais lattice,

*Λ*

^{*}its reciprocal lattice and

*V*the volume of the primitive cell. This property may also be seen by replacing the outgoing free-space Green's function by the incoming one (i.e. its complex conjugate) to obtain (3.2), (3.1) and (3.4). Also, the limiting procedure described in §3

*a*(as an alternative to Ewald's summation) leads to the same conclusion. We were not able to find this fact in the literature except for a particular case in McPhedran

*et al*. (2005).

### (a) A connection with Ewald's method

Ewald (1921) used an integral representation of the free-space Green's function in dimension *d*=3,where *Γ* is a suitably chosen contour in the complex plane so that the integral is well defined. Let us consider a different limiting procedure, similar to that in (2.2), to obtain Ewald's result. Instead of integrating along a contour, we add an imaginary part to *κ*, *κ*+i*λ* with *λ*>*κ* and consider(3.11)with integration over the positive real axis. The expression on the left-hand side of the formula yields the free-space Green's function as *λ*→0, whereas the integral on the right-hand side is well defined only for *λ*>*κ*. However, owing to analytic dependence on *λ*, it is possible to use the integral in (3.11) to obtain an expression for the quasi-periodic Green's function. We proceed with the derivation for dimension *d*=3 and note that in dimension *d*=2 we may follow the same steps but starting (for *λ*>*κ*) withinstead of (3.11). Similar integrals are available in any dimension *d*.

The fact that the integral in (3.11) is well defined for *λ*>*κ* may be seen using the primitiveand, using (Abramowitz & Stegun 1970, eqn (7.1.16)) to evaluate the limits for *λ*>*κ*,andAs a starting point to construct the quasi-periodic Green's function, we use (3.11) and for *λ*>*κ* consider(3.12)As in Ewald (1921), we introduced a real parameter *η*>0 to split the region of integration into two intervals *t*∈(0, *η*) and *t*∈(*η*, ∞).

In the second term in (3.12), we set *λ*=0 since the integral is convergent for all *λ*≥0. Thus, we obtain(3.13)We note that explicit integration yieldsan expression that may be found in some numerical procedures for Ewald's method. We note that this formula requires appropriate modifications in computing contributions of the error function to avoid loss of accuracy. We instead write(3.14)where we changed the integration variable from *t* to *s*, *t*=e^{s}/2. This expression coincides with (3.2) for *d*=3.

In the first term in (3.12), we exchange the order of summation and integration since *λ*>*κ*. We then use the Poisson summation formula in proposition 2.1 to obtainwhere *Λ*^{*} is the reciprocal lattice. By again switching the order of summation and integration, we arrive atDenoting the right-hand side aswe observe that *G*_{Fourier}(** x**,

*λ*) is an analytic function in

*λ*and, thus, may be extended to the region

*λ*>0. This leads us to consider the quasi-periodic Green's function defined by the limiting procedure(3.15)valid for any real parameter

*η*>0. To evaluate the limit, we have (e.g. Gel'fand & Shilov 1964, ch. III, §1.3)For

*κ*>0, we arrive at(3.16)The sum involving the delta function vanishes provided , where and, under this assumption, we obtain (3.1).

## 4. Fast convolutions with Green's function

Representation of the quasi-periodic Green's function as a sum of two rapidly convergent series (3.1) and (3.2) yields a fast and accurate algorithm for its application as a convolution. We truncate these series and obtain a separated representation by approximating the integral in (3.2) via a sum of Gaussians. Using the resulting approximation of Green's function, we prove an accuracy estimate (in operator norm) for its application. We then present the algorithm to apply the operator, and estimate its computational complexity. We illustrate the algorithm by presenting several examples.

### (a) Approximation of Green's function

Let us outline how we obtain an approximation of the quasi-periodic Green's function (3.4).

Owing to the exponential decay of the terms in *G*_{Fourier}, we truncate the Fourier sum(4.1)where we select parameters *η*>0 and *b*>0 so that the contribution of the discarded terms is less than the desired accuracy *ϵ*.

For *G*_{spatial} we perform a similar truncation again using the exponential decay of its terms and, in addition, construct an approximation of *F*_{sing} in (3.3) as a sum of Gaussians. For a fixed parameter *η* and given accuracy *ϵ*, we select *a*>0 to truncate the sum (3.2) asso that the contribution of the discarded terms is less than *ϵ*. Then, for fixed *κ*, we approximate *F*_{sing} as in Beylkin *et al*. (submitted) using a discretization of the integral. Thus, we obtain an approximation of *F*_{sing} as a sum of Gaussians,(4.2)where *σ*_{j}>0 and *q*_{j}>0. The weights *q*_{j} depend on the dimension *d* and the parameter *κ* (see Beylkin *et al*. (submitted) for details). Using (4.2), we approximate *G*_{spatial} as(4.3)Combining (4.1) and (4.3), the quasi-periodic Green's function is approximated as(4.4)We note that there are two sources of error in this approximation: (i) a truncation error due to replacing infinite series by finite sums and (ii) an approximation error introduced by (4.2). Owing to the exponential decay of the terms in both series, the number of significant terms depends only logarithmically on the desired accuracy.

We compute convolutions with in the Fourier domain as(4.5)where(4.6)and convolutions with in the spatial domain as(4.7)We show in dimension *d*=2,3 that this approximation of (3.4) by (4.4) yields accurate convolutions in the operator norm.

*For any ϵ*>0*,* *we may choose the splitting parameter η, the Fourier truncation parameter b in* (*4.1*) *and the spatial truncation parameter a in* (*4.3*) *so that**for* *,* .

We note that the parameters *η*, *b* and *a* are interdependent and their selection is discussed in §4*b*.

Using Minkowski's inequality for convolutions (e.g. Grafakos 2004, p. 20), we haveWe may choose *η*>0 and *b*>1 so thatand, thus,(4.8)We now estimate the spatial error by(4.9)Next, we may choose *a*>0 so that the integrand in the second term satisfiesfor ** x**∈

*D*and, thus,(4.10)In what follows, let us first consider dimension

*d*=3. To estimate the first term in (4.9), as in Beylkin

*et al*. (submitted), we construct the spatial approximation

*S*

_{sing}for accuracy

*ϵ*

_{1}>0 and range parameter 0<

*δ*≤diam(

*D*)/2,(4.11)We use (3.2) and (4.3), where we split the

**=0 term to estimate(4.12)We estimate the first term in (4.12) by considering a ball of radius**

*l**R*=diam(

*D*)/2 circumscribing the primitive cell and using (4.11),(4.13)The second term in (4.12) is estimated asSetting , the minimum diagonal distance in the parallelepiped defined by the lattice vectors, we have(4.14)where

*N*

_{a}is the number of lattice points within the ball . Combining (4.10), (4.13) and (4.14) we have(4.15)For a given

*ϵ*, we select and choose

*ϵ*

_{1}in (4.11) so that . Together with (4.8), we obtain the result for dimension

*d*=3.

For dimension *d*=2, the proof is similar except that we use the spatial approximation (see Beylkin *et al*. submitted)instead of (4.11). ▪

### (b) Choice of the splitting and truncation parameters

The splitting parameter, *η*, controls the rate of decay of the terms in (3.1) and (3.2) and thus, for a given accuracy, determines the number of terms to be retained in each sum. Instead of choosing *η* directly, we choose the Fourier truncation parameter *b*>1 in our approximation (4.1), so that for a given accuracy *ϵ* and *κ*, we set(4.16)With this selection of *η*, note that in (3.1) the discarded terms satisfyWith *η* given by (4.16), we now select the spatial truncation parameter *a* so that the contribution of the discarded terms in (3.2) is below the desired accuracy.

Although we only require *b*>1, in practice the choice of this parameter does depend on *κ* and *ϵ*. For moderate size *κ* we select *b*∼3, for large *κ* we may select a smaller *b* and for small *κ* we need to choose *b* larger.

Different choices of *η* have been made in several papers considering Ewald's summation (e.g. Catti (1978) or Jordan *et al*. (1986) for *κ*=0). We would like to point out (see also Moroz 2006; Oroskar *et al*. 2006 or Beylkin *et al*. submitted) that some choices of *η* may induce numerical cancellation resulting in a loss of accuracy. For example, choosing *η* too small leads to (3.1) and (3.2) to be large simultaneously and to have opposite signs for .

### (c) Algorithm for convolution with the quasi-periodic Green's function

We describe an algorithm and estimate its complexity for computing volumetric convolutions with the quasi-periodic Green's function approximation (4.4). We assume that the input function and its Fourier transform (4.6) are given, and we are free to discretize them as needed. In the description of the algorithm to computewe refer to *f* and *g* as the input and the output function, respectively. We want to compute this convolution for any given accuracy *ϵ*.

Initialization:

*Truncation of the Fourier sum.*For a fixed*κ*and a given accuracy*ϵ*, we select*b*that determines*η*in (4.16).*Truncation and approximation of the spatial sum.*For a fixed*κ*and a given accuracy*ϵ*, we construct*S*_{sing}as an*N*-term Gaussian approximation (4.2). We denote by*N*_{a}the total number of lattice points that satisfy . We note that (see Beylkin*et al*. submitted) and due to the exponential decay of (3.3).*Discretization of the input function.*We use the multiresolution algorithm in Beylkin*et al*. (2008) to adaptively discretize the input function using a tensor product basis with*p*scaling functions per dimension. If*N*_{box}is the total number of boxes used to represent the input function with accuracy*ϵ*, then the total number of input points is . In practical applications, we choose since it improves the overall performance. Thus, we have . We note that it is not hard to construct examples of functions for which an adaptive representation offers no advantage; in such cases, the number of points is due to the required Nyquist sampling rate. Thus, in the worst case, we have .*Initialization of the output function.*The output function, a sum of spatial and Fourier contributions, is evaluated on a user chosen set of*N*_{out}points. While the spatial contribution may retain an adaptive structure if we use the algorithm from Beylkin*et al*. (2008), the Fourier contribution results in (*κ*^{d}) points due to the required Nyquist sampling rate. Thus, unless there are special circumstances, . Again, in the worst case we have .

Applying the operator:

*Convolution with*. Using the algorithm in Beylkin*et al*. (2008), the complexity of applying in (4.7) is . Alternatively, the fast Gauss transform (see Greengard & Strain 1991; Strain 1991; Greengard & Sun 1998) may be used, which results in a similar computational complexity. Although*p*.*N*is formally estimated as , we note that within the range of parameters we experimented with, this product behaves effectively as a constant (the overestimation is, in part, due to the fact that the algorithm in Beylkin*et al*. (2008) does not use all Gaussian terms on all scales). Note that in (4.7) the term=0 dominates the computational cost since this is the only term contributing to fine scales in a multiresolution representation of the operator. With these caveats, the computational complexity of computing (4.7) is , where*l**C*_{3}is a constant.*Convolution with*. We evaluate the Fourier transform of the input function at the reciprocal lattice points within the sphere and denote by*N*_{F}their total number. We note that due to the exponential decay of the terms in (3.1). Given a set of locationsto evaluate (4.5), we use the USFFT (Dutt & Rokhlin 1993; Beylkin 1995; Lee & Greengard 2005) to evaluate the trigonometric sum. Thus, the computational complexity is , or , where*x**C*_{4}is a constant.

We note that the performance of both, the spatial and Fourier, components of our method has been tested and timed independently in the references mentioned above. We would like to add that, in some applications, the semi-analytic nature of our approximation may allow for additional savings.

### (d) Examples

We start by comparing values of the quasi-periodic Green's function computed using in (4.4) and those of its alternative representation in McPhedran *et al*. (2000). The two-dimensional quasi-periodic Green's function is written in McPhedran *et al*. (2000) as(4.17)where and the coefficients are computed as lattice sums. We note that the representation in (4.17) allows us only to evaluate Green's function and does not provide an algorithm for its application as an operator. By contrast, our approach treats Green's function as an operator and constructs an approximation that yields a fast and accurate algorithm for its application. For the purpose of comparison, we implemented the evaluation of Green's function in (4.17) by computing the coefficients in (4.17) as lattice sums, writing . We use (McPhedran *et al*. 2000, eqn (17)) to compute and (Linton 1998, eqns (2.49), (2.53) and (2.54)) to compute .

In figure 1, we display the error between (4.17) and our approximation in (4.4) constructed for accuracy . We note that the discrepancy near *r*=0 is due to our method of approximating *G*_{q} and does not affect its application as an operator (beyond accuracy ) as is demonstrated in proposition 4.1.

Next we verify accuracy of our algorithm by considering the quasi-periodic function(4.18)with parameters *α*=300, ** k**=(1/3, 4/7),

*r*_{1}=(0, 0),

*r*_{2}=(1/10, 1/10) and

*r*_{3}=(−3/20, 1/10), where the sum is over the square lattice generated by the lattice vectors and . We treat

*u*as a solution of (1.4) and (1.5) and analytically compute the corresponding right-hand side in (1.4),(4.19)We construct an approximate two-dimensional quasi-periodic Green's function with

*κ*=30, and apply it to

*f*so that we can compare the result with the exact solution

*u*. The parameter

*α*was chosen so that the Fourier transform of the function in (4.19) remains significant well beyond the disc of radius

*κ*. Such choice allows us to test both the spatial and Fourier parts of the algorithm. In figure 2, we display the absolute error plotted along the diagonal of the primitive cell. Green's function was approximated with

*ϵ*=10

^{−11}, whereas the

*L*

^{2}-norm of the solution is and that of the right-hand side is . This result agrees with the estimate in proposition 4.1.

Next, we illustrate the results of convolving with several quasi-periodic Green's functions. In figure 3, we illustrate the application of a two-dimensional quasi-periodic Green's function to a delta function. The motivation for presenting this example is twofold: (i) to demonstrate that our approach is applicable to functions whose Fourier transforms have slow decay and (ii) to illustrate Green's function itself. In figure 4, we display the result of convolving a periodic Green's function with a fairly complicated function with jump discontinuities. We also display cross sections of the (periodic) output function.

## 5. Green's functions with boundary conditions on simple domains

We now have the necessary tools to construct Green's functions that incorporate boundary conditions on simple domains by extending our results for the quasi-periodic Green's function (3.4). We note that although the resulting integral operators are no longer convolutions, the algorithm for applying these Green's functions is similar to that for the quasi-periodic Green's function. The application of Green's functions satisfying Dirichlet, Neumann or mixed boundary conditions is again split between the spatial and the Fourier domains. In the spatial domain, we use separated representations involving Gaussians and in the Fourier domain apply a simple combination of multiplication operators.

For ease of notation, we consider the two-dimensional case with Dirichlet boundary conditions on the primitive cell . We construct these Green's functions using the periodic Green's function (with 2*κ* instead of *κ*), satisfyingand (1.3) with ** k**=0. We note that the formal description of the periodic Green's function in this case is of the formsince, in (1.6), the sum associated with the imaginary part of the free-space Green's function is zero, and .

We write *G*_{p} via the sum of two rapidly convergent series in (3.4),We obtain Green's function with Dirichlet boundary conditions on *D* as(5.1)For ** x**≠

**, we have since each of the four summands in (5.1) is a Helmholtz Green's function with parameter 2**

*y**κ*. The only singularity is at

**=**

*x***, in which case the first term in (5.1) yieldsSince**

*y**G*

_{p}is periodic with period one and is even, the terms in (5.1) cancel each other on the boundary so that

*G*

^{D}satisfies the Dirichlet boundary conditions, and .

Following the approach in §3, we split (5.1) between the spatial and the Fourier domains and then approximate these components. As in §4*a*, we approximate the spatial part by a sum of Gaussians. For a desired accuracy *ϵ* and fixed *η*, we select *a*>0 to satisfy (4.10) and construct *q*_{j} and *σ*_{j} for *j*=1, …, *N* in (4.11) to obtain the separated representation(5.2)where(5.3)Thus, the application of the operator (5.2) separates along each direction and we computewhich may be accelerated further using fast algorithms described in §4.

In the Fourier domain, for a desired accuracy *ϵ* and fixed *η*, we select *b*>1 to satisfy (4.8) and obtain(5.4)We apply this operator as(5.5)where is given in (4.6). We use USFFT to evaluate (5.5) as in §4*c*.

As described by Keller (1953), the method of images in dimension *d*=2 yields Green's function with prescribed boundary conditions for four bounded regions: (i) rectangle, (ii) equilateral triangle, (iii) isosceles triangles with angles , and (iv) right triangle with angles . As an example of incorporating the Neumann boundary conditions on *D*, we havewhere and , where , *i*=1, 2.

We note that we can mix Dirichlet and Neumann boundary conditions since it requires only appropriate sign changes in the previous construction.

The construction of Green's functions with Dirichlet or Neumann boundary conditions on *D* in dimension *d*=3 is completely analogous to the two-dimensional case and is composed of a combination of eight terms. Importantly, their approximations have the same form in all dimensions. For example, in the spatial domain the approximation of Green's function with Dirichlet boundary conditions is given bywhere *q*_{j} are described in (4.2) and *S*_{j,n} in (5.3). Similarly, we have an analogue of (5.4),which we apply as a multiplication operator in the Fourier domain. In arbitrary dimension *d*, we havewhere are given in (5.3) andWe note that in order to apply Green's function in higher dimensions, we also need to use a separated representation for the input functions (see Beylkin & Mohlenkamp 2005). Green's function with Neumann boundary conditions on *D* has the same form and differs only by changing the sign of appropriate terms. As a result, we may use essentially the same algorithm to apply these operators. To summarize, the results of this section yield fast adaptive solvers for the Helmholtz equation for a variety of boundary conditions.

## 6. Conclusion and remarks

In this paper, we extend the approach in Beylkin *et al*. (submitted) for the free-space Helmholtz Green's function to approximate and apply Green's functions, which incorporate quasi-periodic Dirichlet or Neumann boundary conditions. The key features of these fast algorithms are: (i) the splitting of application of operators between the spatial and the Fourier domains, (ii) the use of separated representations, and (iii) the ability to achieve a finite, arbitrary accuracy. Algorithms with the last two features have been developed for non-oscillatory kernels and have been used to solve problems in quantum chemistry (see Harrison *et al*. 2003, 2004; Yanai *et al*. 2004*a*,*b*). Since these algorithms for oscillatory and non-oscillatory kernels may be considered within the same framework, we intend to build a unified software framework for their application. We expect further development in this direction. In all cases, we obtain representations of Green's functions that lead to fast adaptive solvers for corresponding problems.

Our approach (with minor modifications) is also applicable to the case *κ*=0. However, using multiresolution, both the interpretation and the application of the operator may be kept entirely in the spatial domain and we plan to consider this case separately.

A natural application of the quasi-periodic Green's function is in the computation of band gaps in crystal structures. We plan to investigate these applications with particular attention to potentials (indices of refraction) with singularities (discontinuities) since, in such cases, the efficiency of our algorithms does not degrade significantly.

We note that our method extends to problems where the lattice dimension is less than the dimension of the embedding space (sometimes referred to as gratings), which will be described elsewhere.

Finally, we note that our results shed new light on Ewald's approach of splitting between spatial and Fourier domains, which we use as a tool to obtain semi-analytic, separated representations for Green's functions.

## Acknowledgments

This research was partially supported by NSF grant DMS-0612358, DOE/ORNL grant 4000038129, DOE grant DE-FG02-03ER25583 and AFOSR grant FA9550-07-1-0135.

## Footnotes

- Received April 18, 2008.
- Accepted July 22, 2008.

- © 2008 The Royal Society