Method of lines transpose: an efficient A-stable solver for wave propagation

Building upon recent results obtained in [7,8,9], we describe an efficient second order, A-stable scheme for solving the wave equation, based on the method of lines transpose (MOL$^T$), and the resulting semi-discrete (i.e. continuous in space) boundary value problem. In [7], A-stable schemes of high order were derived, and in [9] a high order, fast $\mathcal{O}(N)$ spatial solver was derived, which is matrix-free and is based on dimensional-splitting. In this work, are interested in building a wave solver, and our main concern is the development of boundary conditions. We demonstrate all desired boundary conditions for a wave solver, including outflow boundary conditions, in 1D and 2D. The scheme works in a logically Cartesian fashion, and the boundary points are embedded into the regular mesh, without incurring stability restrictions, so that boundary conditions are imposed without any reduction in the order of accuracy. We demonstrate how the embedded boundary approach works in the cases of Dirichlet and Neumann boundary conditions. Further, we develop outflow and periodic boundary conditions for the MOL$^T$ formulation. Our solver is designed to couple with particle codes, and so special attention is also paid to the implementation of point sources, and soft sources which can be used to launch waves into waveguides.


Introduction
Numerical solutions of the wave equation have been an area of investigation for many decades. The wave equation is ubiquitous in the physical world, arising in acoustics, electromagnetics, and fluid dynamics. Our main interest is in electromagnetic wave propagation, where traditional finite difference methods such as the finite-difference time-domain (FDTD) method are often used to solve Maxwell's equations. When the classical Yee scheme is used, the Courant-Friedrichs-Lewy (CFL) stability criterion restricts the time step to scale with the smallest cells in the domain. This becomes computationally prohibitive in a variety of interesting problems, such as electromagnetic scattering or waveguide design, where complex geometries need to be embedded in a Cartesian mesh, leading to small spatial cells near the boundaries. Alternatively, problems with multiple and disparate temporal or spatial scales, such as those presented by plasma simulations, require time steps and mesh spacings which are on the order of the shortest temporal and spatial scales. Due to the large number of charged particles in a typical simulation [4], and the high dimensionality of their distribution space, time steps are computationally prohibitive, and the amount of them must be minimized.
As a result, A-stable Maxwell solvers have been developed, which remove the CFL restriction and restore the ability of the user to define a time step which is based on the physical problem, rather than on its spatial discretization. Notable advances include the introduction of the FDTD-ADI [14,15,23,24] algorithm, as well as several semi-implicit time split schemes. However, it remains difficult to preserve both accuracy, and A-stability when non-rectangular domains are required.
Alternatively, Maxwell's equations can be reformulated (e.g., using the scalar and vector potential formulation), so that each field component independently satisfies the second order wave equation. Method of lines (MOL) formulations of the wave equation are well-studied, and often lead to conditionally stable schemes. But when the discretization is first performed in time, following the method of lines transpose (MOL T ), the wave equation is found at discrete time levels by solving a semi-discrete boundary value problem. In [1,2], the MOL T is used to produce an exact integral formulation of the wave equation, where the solution is determined by a convolution over the domain of dependence against a space-time Green's function G d (x, t), in dimensions d = 1, 2, 3. However, more commonly the semi-discrete solution is obtained, in which the modified Helmholtz equation must be solved, and the corresponding semi-discrete Green's function G(x, ∆t) is the Yukawa, or modified Helmholtz kernel. The resulting boundary integrals methods can then be solved at each time step using fast summation algorithms, such as tree-codes [3,12,20,21], or the fast multipole method (FMM) [17,10,13,19,16,20,11]. These methods often scale as O(N ) or O(N log N ), but require substantial storage or precomputing stages.
It is worth noting however that in one spatial dimension, the modified Helmholtz kernel is a simple decaying exponential. This fact has been combined with alternate dimension implicit (ADI) splitting to achieve an A-stable wave solver with computational complexity of O(N log N ) [6,22], and O(N ) [8,7], respectively. The scaling presented in [6,22] is due to the Fourier continuation method, which makes use of the fast Fourier transform to compute a periodic extension of the the convolution integral, where the period is taken sufficiently longer than the domain, so that boundary conditions can be satisfied.
In our approach [7,8,9], the one-dimensional convolution is instead computed strictly over the computational domain, using polynomial interpolation. In [8], the solution is proven to be A-stable, second order accurate, and that the matrix formed by the discrete convolution can be applied in O(N ) operations. However, the discrete convolution matrix was formed using piecewise linear integration, and it was found that a Lax-type correction was required to ensure convergence of the scheme in the semi-discrete limit ∆t → 0, where ∆x is fixed. This issue was addressed in [9], where spatial discretization was extended to orders M ≥ 2, on non-uniform grids while retaining O(N ) complexity. In [7], the scheme was extended to higher orders in time using a novel approach, successive convolution. These higher order schemes were proven to be A-stable.
The purpose of this present work is to develop a robust embedded boundary approach for complex geometry for the MOL T formulation. Because this paper is focused on developing a range of boundary conditions for the wave equation, we limit our attention to the 2nd order accurate solver. In addition to addressing standard closed (i.e., Dirichlet, Neumann and periodic) boundary conditions, we also develop an open, or outflow boundary condition which is suitable for our implicit solver. The method is extended to higher spatial dimensions using the factorization developed in [7], which is similar to but different from the traditional ADI splitting. This means that solutions are computed line-by-line along dimensional sweeps leveraging the O(N ) 1D solver to construct a high dimensional O(N ) implicit method.
The rest of this paper is laid out as follows. In section 2, we use the method of lines transpose (MOL T ) to reformulate the 2D wave equation as an implicit, semi-discrete modified Helmholtz equation. The Helmholtz operator is inverted analytically using dimensional splitting, and we recast the solution as a series of one-dimensional boundary integral equations. We specifically show how to obtain a fully discretized, second order accurate solution by performing spatial quadrature, using a fast O(N ) convolution algorithm.
In section 3, we modify the time-centered scheme of [7] to introduce artificial dissipation. This will be used to stabilize the embedded boundary method in our 2D algorithm, presented in section 4. We demonstrate our solution to be fast, second order accurate, and A-stable, even on non-rectangular domains with numerical results presented in section 5. We conclude the body of the paper with several remarks in section 6.
We also include in Appendix A, a table summarizing the main algorithmic aspects of our wave solver. In Appendix B, we show how point sources, which may represent charged plasma particles, or soft point sources launched into a waveguide, can be implemented.

Integral solution using MOL T
We now develop a dimensionally split algorithm for solving the initial value problem with consistent boundary conditions. We utilize the method of lines transpose (MOL T ) to perform a second order accurate temporal discretization, as was shown in [7]. In particular, the discretization is time-centered, and implicit in the spatial derivatives and β ≤ 2 is chosen to enforce A-stability [7]. Inversion of this operator leads to a boundary integral formulation, where the Green's function is the Yukawa potential, defined for 2D in terms of the modified Bessel function K 0 (r). However we will employ our previously developed [7,8] splitting algorithm, which has O(N ) complexity Here the subscripts denote the spatial component of univariate modified Helmholtz operators, The modified Helmholtz operators are formally inverted using the Green's function. We define convolution with this Green's function by the integral operator so that where the coefficients A and B of the homogeneous solution are determined by applying boundary conditions. A similar definition holds for L −1 y . Formally inverting both operators to the right hand side we find the explicit equation where the multidimensional operator is now . Since the application of L −1 x is done for fixed y, and vice versa for L −1 y , the operator D xy can be constructed in a line-by-line fashion, similar to ADI algorithms. It was also proven in [7] that the scheme (6) is A-stable, for 0 < β ≤ 2. Remark 1. In the continuous case, L x and L y (and their inverses) commute, so D xy = D yx . However in practice the discretize operators will not commute, and some small spatial anisotropy is introduced. This can be controlled by applying both operators D xy and D yx , and then averaging the result.
Remark 2. Similar to a traditional ADI formulation of the wave equation, this factorization produces a fourth order splitting error term, which can be compensated for by adding a term to the right hand side Note our use of the identity (L x − 1)(L y − 1) = ∂ xx ∂ yy /α 4 . Formally inverting both operators to the right hand side, the solution can be rearranged and found as where the convolution operator C is the operator defined in [7] as This form was used in [7] to achieve schemes of higher order through successive convolution, where removing splitting errors is of paramount concern.
2.1. Fast convolution algorithm for the integral solution. In [8], the particular solution (4) was discretized in space using the weighted midpoint and trapezoidal rules, which amounted to replacing u with a piecewise constant and linear approximation, respectively. More recently [9], we have detailed the spatial discretization of I x to arbitrary order, while providing a means for its rapid evaluation. In this work we focus on the second order accurate implementation of this algorithm, and reiterate the relevant details. This particular solution is first decomposed into a left and right oriented integral, split at y = x so that where These quantities can be updated locally, using exponential recursion The recursive updates (11) and (12) are exact (in space), and making δ L and δ R small (typically, δ L = δ R = ∆x) effectively localizes the contribution of the integrals. Based on these observations, we now outline the fast convolution algorithm. In [9], we derived a spatial quadrature of general order M ≥ 2 for irregular grids. Here, we will utilize a second-order accurate method on a uniform grid, which we explicitly develop.
Consider the domain (a, b) discretized by uniform grid points x 1 = a < x 2 < · · · < x N +1 = b of width ∆x = b−a N = x j+1 − x j , j = 1, ..., N . Suppose, given a function f compactly supported in (a, b), we are to evaluate the convolution operator at each grid points, that is compute We proceed by evaluating the local integrals J L j = J L [f ](x j ) and J L j = J L [f ](x j ), as in (11) and (12). The local integrals may be evaluated with quadrature, or, if possible, analytically. A second-order accurate quadrature is given by where defining ν = α∆x and d = e −ν , the quadrature weights are given by We summarize our fast method in Algorithm 1.

Artificial dissipation
As will be discussed in section 4.2.2, it is necessary to include some artificial dissipation in the numerical scheme to maintain stability with embedded boundary methods for Neumann boundary conditions. We first present a version of the wave solver based on a backwards difference formula (BDF) time discretization, leading to what we call a diffusive scheme, which is dissipative. This method has a larger truncation error than a centered scheme, does not possess a means to tune the level of dissipation, and also has an implicit source term (at time level n + 1), which is problematic for application in the context of particlein-cell (PIC) methods for the simulation of plasmas. The second-order centered (dispersive) scheme given above is therefore preferable, but are non-dissipative in their original forms. We give a method for adding tunable artificial dissipation terms into the centered scheme, while maintaining A-stability.
3.1. Diffusive wave solver. We substitute the following backward difference formula (BDF) discretization: into the wave equation 1 c 2 u tt − ∇ 2 u = S(x, t). Rearranging, defining α = √ 2/(c∆t) and dividing by α 2 gives the semi-discrete scheme This method is A-stable and dissipative, but does not possess a mechanism for tuning the dissipation, has an inconvenient implicit source term (at time level n + 1), and typically has a larger truncation error compared to the centered scheme.

Artificial dissipation in centered schemes.
3.2.1. Artificial Dissipation in 1D. We give a modified form of the centered version of the implicit wave solver with tunable artificial dissipation that retains the property of unconditional stability. We let denote a small artificial dissipation parameter, and D Ignoring sources, we have the second order scheme with dissipation, We now prove the unconditional stability of this scheme for prescribed values of β. As in [7], we pass to the high-frequency limit.
We obtain the Von Neumann polynomial ρ 2 − (2 − β 2 )ρ + (1 − ). We can check that the roots of this polynomial will be complex if 0 < β ≤ 2 + 2 √ 1 − , and that in this case the roots satisfy which shows both the stability and dissipative nature of the scheme. We note that the maximum allowed value of β is slightly smaller than what is allowed by the corresponding scheme without dissipation. A more detailed analysis shows that the effective damping rate is k 2 k 2 +α 2 2 , meaning that high frequencies are more rapidly damped than low frequencies.
3.2.2. Artificial Dissipation in 2D. Using the notation defined in [7], and again ignoring sources, we have the second order scheme with dissipation, For further details on these operators, see [7]. Numerical experiments indicate that the 2D scheme with artificial dissipation is indeed unconditionally stable with the same maximum value of β as with the 1D schemes.

Boundary conditions
We will now discuss the implementation of boundary conditions. Since our algorithm is dimensionally split, we first develop the boundary conditions in one spatial dimension where the solution (6) reduces to and where we consider the following boundary conditions Periodic Once these boundary conditions have been derived, we use them to build a boundary solver in 2D.

4.1.
Boundary conditions in one dimension. In 1D, this homogeneous solution requires the determination of two coefficients from the imposed boundary conditions and the endpoint values of the particular (integral) solution by solving a 2 × 2 linear system. We now show how to impose several common boundary conditions in 1D. These methods are extended to the 2D case in Section 4.2.
4.1.1. 1D Dirichlet boundary conditions. Let us begin with Dirichlet boundary conditions (25). Evaluating the semi-discrete solution (24) at x = a and b, we find which, after solving for the unknown coefficients can be written as . Homogeneous boundary conditions are recovered upon setting U L (t) = U R (t) = 0. Solving the resulting linear system for the unknowns A n and B n gives 4.1.2. 1D Neumann boundary conditions. For Neumann conditions, first observe that all dependence on x in the integral solution (5) is on the Green's function, which is a simple exponential function. Using this, we obtain the following identities Now, differentiating the semi-discrete solution (24), and applying the Neumann boundary conditions (26) at x = a and b yields which, after solving for the unknown coefficients can be written as Upon solving the linear system we obtain Remark 3. The cases of applying mixed boundary conditions at x = a and b are not considered here, but the details follow from an analogous procedure to that demonstrated above.

4.1.3.
1D Periodic boundary conditions. We impose periodic boundary conditions, by assuming that u n (b) = u n (a), u n x (a) = u n x (b), n ≥ 0. Enforcing this in the semi-discrete solution (24) then yields where we have used the identity (30) applied to derivatives of I. Solving this linear system is accomplished quickly by dividing the second equation by α, and either adding or subtracting it from the first equation, to produce 4.1.4. 1D Outflow boundary conditions. When computing wave phenomena, whether we are interested in finite or infinite domains, it is often the case that we must restrict our attention to some smaller subdomain Ω of the problem, which does not include the physical boundaries. We say that Ω is the computational domain, and that the boundary ∂Ω is the non-physical, or artificial boundary. Under these circumstances, it is necessary to enforce an outflow, or non-reflecting boundary condition, which allows the wave to leave the computational domain, but not incur (non-physical) reflections at the artificial boundary. For this reason, let us consider the free space solution but where we are only interested in evaluating this expression for x ∈ Ω = [a, b]. Then the contributions can be decomposed as where the homogeneous coefficients are which we observe do not depend on x. Since the coefficients A an B are the contributions of the integral to the left and right of [a, b] respectively, they can be thought of as transmission conditions (rather than boundary conditions). We make use of this fact to develop outflow boundary conditions, and it will serve as a key idea in our planed follow-up work on a domain decomposition algorithm and multi-core computing with our implicit wave solver. For the one-dimensional wave equation the exact outflow boundary conditions (28) turn out to be local in space and time. We emphasize that this is only the case in one spatial dimension, but we shall utilize this fact to obtain an outflow boundary integral solution from the integral equation (33). We extend the support of our function to (−∞, ∞), and extend the definition of the outflow boundary conditions to the domains exterior to [a, b] Next, assume the initial conditions have some compact support; for simplicity we will take this support to be Ω 0 = [a, b]. Then after a time t = t n , the domain of dependence of u n (x) is Ω t = [a − ct n , b + ct n ], since the propagation speed is c. Now the free space solution (33) becomes with coefficients At first glance, these coefficients are not at all helpful, as they require computing integrals along spatial domains which not only are outside of the computational domain, but also grow linearly in time. However, we will now make use of the extended boundary conditions to turn these spatial integrals into time integrals, which exist at precisely the endpoints x = a and b respectively. Consider first x > b. By assumption, this region contains only right traveling waves, u(x, t) = u(x − ct), and by tracing backward along a characteristic ray we find Thus, and so B n is equivalently represented by a convolution in time, rather than space. Now, knowing the history of u at x = b is sufficient to impose outflow boundary conditions. Furthermore, we find in analog to equation (12), a temporal recurrence relation due to the exponential where β = αc∆t, by definition (3). Thus, the coefficient B n , which imposes an outflow boundary condition at x = b, can be computed locally in both time and space. To maintain second order accuracy, we fit u with a quadratic interpolant and integrate the expression analytically to arrive at where In this outflow update equation (41), the quantities u n+1 (b) and B n are both unknown. In order to determine these values, we must also evaluate the update equation for u n+1 (24) at We now use these two equations to solve for u n+1 (b), and eliminate it from the outflow update equation (41), so that where Remark 4. While this procedure could be avoided by omitting u n+1 (b) in the interpolation stencil, it turns out to be necessary to obtain convergent outflow boundary conditions.
Likewise, upon considering x < a, we find Solving the resulting linear system produces where w Out Boundary conditions in two dimensions. We now describe our approach for imposing boundary conditions up to second-order accuracy in two dimensions for our MOL T formulation of the wave equation. Boundary conditions must be supplied for the intermediate sweep variable w. Since w = u + O(c 2 ∆t 2 ), for second order accuracy it suffices for w to inherit the boundary condition imposed on the main solution variable u. In the case of rectangular, grid-aligned boundaries, the 1D boundary correction terms can be imposed   Thus, a line y = y k , which is discretized in x, and has endpoints a k and b k , determined as the intersection of the line with the boundary curve. Analogously we define endpoints of the line x = x j as c j and d j , and now have for 1 ≤ k ≤ N y and 1 ≤ j ≤ N x respectively. If the boundary is defined using a level set function C(x, y) = 0, then the points x = a k , b k can be found by solving C(x, y k ) = 0, and following from the analogous approach, the endpoints y = c j , d j corresponding to x = x j are found.
Once the domain has been discretized and the lines have been defined, it remains to compute the discrete form of the scheme (6). Since the one-dimensional convolution algorithm presented in Section 2.1 is formulated for non-uniform grid points, the embedded boundary points do not affect the implementation of the horizontal (45) and vertical (46) sweeps. Thus, the only point of consideration that remains is that of boundary conditions.
First we consider an intermediate variable w (1) and which, according to the implicit scheme (2), can be seen to satisfy Boundary conditions must be applied to w (1) at the boundary points which terminate each line y = y k , defined as (a k , y k ) and (b k , y k ). Since u n+1 , u n , and u n−1 , will be prescribed at these boundary points, we see that the boundary condition will be of the form 1 − ∂ yy α 2 u n+1 − 2u n + u n−1 + β 2 u n (x, y), and the order of the limit and the partial derivatives commute only when partial derivatives of the boundary data can be constructed. But this is the case only when the tangent line at the boundary is in the y-direction, which is precisely why we must restrict our attention to the aforementioned cases for boundary conditions. Upon introducing a second intermediate variable w (2) (x, y), defined by we make the observation that L y [w (2) ] = w (1) . Boundary conditions are now applied to w (2) along the lines x = x j , at the boundary points (x j , c j ) and (x j , d j ), where now no difficulties remain in considering since the right hand side will be fully prescribed. We summarize this procedure in Algorithm 2.

Algorithm 2 Application of Boundary Conditions in 2 dimensions
(1) Initialize temporary variables w (1) and w (2) , which are the same size u n .
(2) For each horizontal line y = y k , for 1 ≤ k ≤ N y , create the temporary variable w k ](x) = β 2 u n + 1 α 2 S n (x, y k ). The homogeneous coefficients are determined by the fact that which is applied at x = a k , b k .
k (x). The homogeneous coefficients are determined by the fact that which is applied at y = c j , d j .
(5) The dimensional splitting error is corrected by adding

Neumann boundary conditions in two dimensions.
In implementing Neumann boundary conditions for boundary geometries conforming to grid lines, such as a rectangular domain, we can directly impose a two-point boundary correction. One way to extend this method to a general polygonal domain would be to use multiple overset grids, each aligned with a boundary segment, which communicate with the interior grid through interpolation on a ghost cell region, though we do not pursue that approach in this work.
For curved boundaries, an alternative approach is that of an embedded boundary method, which involves determining the Dirichlet values at the endpoints of each x-and y-sweep lines that result in the approximate satisfaction of the Neumann boundary condition (in effect, constructing an approximate Neumann-to-Dirichlet map). We present the implementation of an embedded boundary method for Neumann boundary conditions for the implicit wave solver on a curved boundary geometry. The approach taken here follows the work in [18], which proposes an embedded boundary method for Neumann boundary conditions with a finite difference method for the wave equation. The analysis in that work suggests that, on the continuous level, the modified equations and boundary conditions resulting from typical truncation error terms possess unstable boundary layer solutions, so that the addition of a dissipative term is necessary to achieve a stable method. This is consistent with our experience in the implementation described here, with the embedded boundary method becoming unstable when applied to the non-dissipative dispersive solver, but remaining stable for the diffusive solver, which is dissipative.
In the following, we briefly describe the two-point boundary correction method for a gridaligned rectangular boundary. We then describe the embedded boundary method for a 1D problem. This method requires an iterative procedure, which we show, in the setting of the 1D problem, to be a convergent contraction mapping with a rate of convergence that depends on the CFL number. We describe the implementation of the embedded boundary method in 2D, and finally give numerical results.

Description of the Two-Point Boundary Correction Method.
In a rectangular domain where the boundaries conform to grid lines, it is straightforward to impose the two-point boundary correction terms in a line-by-line fashion, since in this case, the grid lines are not coupled through the normal derivative. As this is a simple extension of the 1D boundary correction algorithm, we do not elaborate further.

Description of the Embedded Boundary Method and Proof of Convergence of the Iterative Solution in 1D.
We consider the situation of a one-dimensional domain {x B < x} with a single boundary point not aligned with the grid points, as displayed in Figure 2. We have grid points x 0 , x 1 , ... with uniform grid spacing x i+1 − x i = ∆x, boundary location x B , and ghost point location x G = x 0 . We define interior points to be any grid points lying within the domain (including the boundary), and exterior points to be any grid points lying outside of the domain. We define a ghost point to be any exterior point for which at least one of the neighboring points x i±1 is an interior point. We neglect the right boundary in the present analysis for simplicity, though it can be extended to the case with both boundaries. We consider applying the diffusive version of the wave solver, having calculated the convolution integral I(x), and now needing to find the value of the coefficient A such that the solution u(x) at the next time step is given by Given the value of the convolution integral and the solution at the ghost point, I G = I(x G ) and u G = u(x G ), respectively, the coefficient may be computed as A = u G − I G . We now describe the procedure for determining the value of the solution at the ghost point, u G , that leads to a solution consistent with homogeneous Neumann boundary conditions to second-order accuracy. We construct a quadratic interpolant using the boundary condition and interior interpolation points x I = x B + ∆s I and x II = x B + 2∆s I , lying in between grid points x m and x m+1 , and x n and x n+1 , respectively. The interpolation distance will be chosen such that ∆x < ∆s I < (3/2)∆x. We define the distances ξ G = x B − x G , ξ I = x I − x B = ∆s I , and ξ II = x II − x B = 2∆s I , and construct a quadratic Hermite-Birkhoff [5] interpolant P (ξ) by imposing the conditions P (0) = 0, P (ξ I ) = u I and P (ξ II ) = u II . We then obtain the following second-order approximation to the ghost point value, given by As the coefficients γ I = , we only need supply second-order accurate approximations to u I and u II to maintain overall second-order accuracy. (The coefficients would be O(1/∆x) in the case of nonhomogeneous Neumann boundary conditions, which would require third-order accurate approximations to u I and u II . For simplicity, we consider only the case of homogeneous Neumann boundary conditions in the present work.) Such approximations may be obtained through linear interpolation, giving , and u j = u(x j ) are the values of the function at the uniform gridpoints for j = m, m + 1, n, n + 1.
x 0 x G x B x 1 x m x m+1 x I x n x n+1 x II Figure 2. Boundary geometry in 1D.
Hence, to determine the ghost point value u G that leads to a solution consistent with homogeneous Neumann boundary conditions, we should solve the following system of equations.
with γ I , γ II , σ I and σ II defined as above, and where I j = I(x j ) is the convolution integral evaluated at uniform grid points for j = m, m + 1, n, n + 1. This system can be solved formally for u G , giving To show that this solution formula is well-defined, we argue that for the relevant values of m, n and ξ G , ξ I , and ξ II . We define d = e −α∆x , and noting that 0 < d < 1, m < n, ξ G < ξ I = ∆s I < ξ II = 2∆s I , γ I > 0, γ II < 0, 0 ≤ σ I ≤ 1, and 0 ≤ σ II ≤ 1, we obtain Now, since ∆x < ∆s I < (3/2)∆x, we can see that it is the case that either m = 1 and n = 2 or 3, or that m = 2 and n = 3. It is then a matter of some simple calculus to check that that the functions f m,n (x) = (4x m − x n+1 )/3 satisfy f m,n (x) < 1 for 0 < x < 1 and the given combinations of m and n. This proves that K < 1 for the relevant values of the parameters, so that the solution for u G is well-defined. We note, however, that K approaches 1 as d approaches 1, that is, as the CFL number becomes large. Thus, we may expect an ill-conditioned system when the CFL number is very large.
To obtain the lower bound on K, we observe In the two-dimensional case, the line-by-line solution method couples the ghost point values, and a general explicit solution formula is impossible to write down. In principle, one may write out and directly solve a linear system to obtain the ghost point values. Instead, we propose an iterative solution method that avoids the formation of a matrix. We now describe this iterative solution method and prove its convergence in the context of the one-dimensional problem described above.
Suppose we have the convolution integral evaluated at the gridpoints, I j , and a k-th iterate for the ghost point value, u k G . Then we may obtain the next iterate by the formulas u k+1 ) where quantities are defined as above. Now, to prove the convergence of the interation, we show it is contractive. Taking the difference of two iterates, we have Hence, the Contraction Mapping Theorem implies that the iteration converges to a unique fixed point (which is the solution given in the formula above). We note again that K approaches 1 as the CFL number becomes large, so that the rate of convergence will become slower for larger CFL numbers.

Description of the Method in 2D.
We now describe the implementation of the embedded Neumann boundary condition in the 2D case. We consider the situation displayed in Figure 3, in which we need to determine the value of our unknown u G at the ghost point location (x G , y G ). In the 2D case, we define a ghost point to be any exterior point (x i , y j ) for which at least one of the neighboring points (x i±1 , y j ) or (x i , y j±1 ) is an interior point. Similarly to the 1D case, we will construct a quadratic Hermite-Birkhoff boundary interpolant P (ξ) along the direction normal to the boundary, which intersects the boundary curve Γ at location (x B , y B ), and supply the interior interpolation point values u I and u II , at points (x I , y I ) and (x II , y II ), respectively, by further interpolation from interior grid points. These points are selected along the normal, in analogy to the 1D case, such that ξ I = |(x I , y I ) − (x B , y B )| = ∆s I and ξ II = |(x II , y II ) − (x B , y B )| = 2∆s I , where we will typically take ∆s I = √ 2∆x. We construct a quadratic Hermite-Birkhoff interpolant P (ξ) by imposing the conditions P (0) = 0, P (ξ I ) = u I and P (ξ II ) = u II . Defining further the distance from the boundary to the ghost point ξ G = |(x G , y G ) − (x B , y B )|, we obtain, as in the 1D case, the following second-order approximation to the ghost point value, given by where the coefficients are defined as γ I = In the 2D case, we find approximations to u I and u II through bilinear interpolation. This is in contrast to [18], who find the intersection of the normal with grid lines, then interpolate along the grid lines. We have also implemented a second-order accurate version of this approach and compared to the bilinear interpolation scheme proposed here. We have found that the two schemes behave similarly, however the bilinear interpolation scheme is slightly more accurate and simpler to code, not requiring to handle separate cases of intersection with horizontal and vertical grid lines. The bilinear interpolation scheme is standard, but we give it here for completeness. If the interpolation point u I lies in a cell with corners (x i , y j ), (x i+1 , y j ), (x i+1 , y j+1 ) and (x i , y j+1 ), then we have the following approximation for u I : . With this interpolation scheme established, we now outline the algorithm for the 2D dimensionallysplit wave solver.
Γ (x G , y G ) (x I , y I ) (x II , y II ) (x B , y B ) Figure 3. Boundary geometry in 2D.
The above interpolation procedure applies regardless of the variety of the wave solver that it is used with, provided that the wave solver has sufficient dissipation to maintain stability. We now describe the rest of the embedded boundary algorithm in the context of the diffusive wave solver, though it it may be similarly applied to the the dispersive scheme with artificial dissipation described above. In analogy to the iteration presented in the 1D case, the 2D iterative algorithm proceeds by using the interpolation scheme to provide values at the ghost points, which in turn provide new values for the boundary correction coefficients, which are then used to update the values at the interior grid points, comprising one full iteration. It should be noted that not all interior grid points need be updated in the iteration, only those near the boundary that lie within the boundary interpolation stencils.
Using values from previous time steps, the initial guess for the interior grid points in the boundary interpolation stencils may be given by extrapolation in time, as u n+1,0 = 2u n −u n−1 (linear extrapolation) or u n+1,0 = 3u n − 3u n−1 + u n−2 (quadratic extrapolation). Either extrapolated initial guess provides a modest reduction in the number of iterations required versus a zero initial guess, with only a slight further reduction in the number of iterations going from linear to quadratic extrapolation. An effective stopping criterion for iteration is |u n+1,l+1 − u n+1,l | ∞ < , where is some chosen tolerance, which may be chosen to be quite small, as the iteration is a fixed point interation. In the numerical example, we choose a tolerance of 10 −15 , and we achieve convergence in less than 40 iterations at a CFL number of 2.
In applying the diffusive version of the wave solver, we assume we have previous time steps u n , u n−1 and u n−2 . We have to solve the modified Helmholtz equation with homogeneous Neumann boundary conditions, c∆t . We apply dimensional splitting to find so we define w = 1 − 1 α 2 ∂ yy u, and noting that w = u + O (c∆t) 2 so that ∂w ∂n = ∂u ∂n + O (c∆t) 2 , we obtain the following approximate system We now suppose our domain is embedded in a uniform Cartesian grid, with horizontal grid lines corresponding to y = y k , 1 ≤ k ≤ N y and vertical grid lines corresponding to x = x j , 1 ≤ j ≤ N x . The embedded boundary algorithm will be applied when calculating the intermediate variable w n+1 in horizontal line sweeps as well as the solution variable u n+1 in vertical line sweeps. The iterations for these two variables are separate; first, the iterative procedure is applied to w to convergence, and then this value of w is used to compute u, and the iterative procedure is applied to u to convergence. However, in each iteration, the grid lines are coupled through the interpolation scheme, so that all grid lines must be iterated together. The overall iterative algorithm is described in 3, with details specified for the iteration on w. The iteration on u is very similar, and so we omit the details.

Algorithm 3 Application of Neumann Boundary Conditions in 2 dimensions
(1) (Initialization of ghost points) Perform the interpolation scheme described above to obtain the values of u n , u n−1 , and u n−2 at the ghost points, which are the endpoints of the horizontal and vertical grid lines.
(2) (Evaluation of particular solution) For each horizontal line y = y k , for 1 ≤ k ≤ N y , with ghost (end) points x = a k and b k find the particular solution w n+1 p,k for the intermediate variable w n+1 k (x) by evaluating the discrete convolution operator for each grid point x j in the horizontal line, including the ghost points. to calculate the particular solution for u n+1 , then apply the bounday correction interation.

Double Circle Cavity.
In this example, we solve the wave equation with homogeneous Dirichlet boundary conditions on a 2D domain Ω which is, as in Figure 4, the union of two overlapping disks, with centers P 1 = (−γ, 0) and P 2 = (γ, 0), respectively, and each with radius R: where | (x, y) | = x 2 + y 2 is the usual Euclidean vector norm, and γ < R.
This geometry is of interest due to, for example, its similarity to that of the radio frequency (RF) cavities used in the design of linear particle accelerators, and presents numerical difficulties due to the curvature of, and presence of corners in, the boundary. Our method avoids the staircase approximation used in typical finite difference methods to handle curved boundaries, which reduces accuracy to first order and may introduce spurious numerical diffraction. As initial conditions, we choose |(x,y)−P 2 | 0.8γ 2 | (x, y) − P 2 | < 0.8γ 0 otherwise and u t (x, y, 0) = 0 for (x, y) ∈ Ω. Selected snapshots of the evolution are given in Figure 5, and the results of a refinement study are given in Table 1. The discrete L 2 error was computed against a well-refined numerical reference solution (∆x = 2.1875 × 10 −4 ); the error displayed in the table is the maximum over time steps with t ∈ [0.28, 0.29]. For this example, R = 0.3, γ = 0.2, c = 1, and the CFL is 2. In both cases, the solution converges to second order. An overlay of the two are shown in Figure 6, demonstrating the close agreement.  Figure 6. Two separate numerical constructions of a Bessel mode are superimposed, demonstrating that the solution on the quarter circle using Neumann boundary conditions is equivalent to that of the full circle.

Bessel Mode with Neumann Boundary Conditions.
Here we present a numerical example of the embedded boundary method for homogeneous Neumann boundary conditions given in Section 4.2.2. We apply the method to a circular domain, for which analytical solutions exist. We consider a radially-symmetric Bessel mode with homogeneous Neumann boundary conditions, with an analytic solution given by where J 0 is the Bessel function of the first kind of order 0, r = x 2 + y 2 , R is the radius of the domain, and Z 0 ≈ 3.8317 is the smallest nonzero root of J 0 (so that ∂u ∂n (R, t) = ∂u ∂r (R, t) = Z 0 R J 0 (Z 0 ) cos Z 0 ct R = 0). In this example, we take radius R = π/2 and wave speed c = 1. An example of the embedded boundary grid used is given in Figure 7. We perform a refinement study with a fixed CFL number of 2, with the results in Figure 8 indicating the expected second-order convergence. We set the iteration tolerance to 10 −15 , and we see convergence of the boundary correction iteration in fewer than 40 iterations. We note some oscillation of the L ∞ error about the line giving second-order accuracy, which we believe to be due to the grid points moving with respect to the boundary through the refinement, causing some variation in the maximum error. 5.4. Periodic Slit Diffraction Grating. In this example, we apply our method to model an infinite, periodic diffraction grating under an incident plane wave. Diffraction gratings are periodic structures used in optics to separate different wavelengths of light, much like a prism. The high resolution that can be achieved with diffraction gratings makes them useful in spectroscopy, for example, in the determination of atomic and molecular spectra. Our numerical experiment, depicted in Figure 9, demonstrates the use of our method with multiple boundary conditions and nontrivial geometry in a single simulation to capture complex wave phenomena.
An idealized slit diffraction grating consists of a reflecting screen of vanishing thickness, with open slits of aperture width a, spaced distance d apart, measured from the end of one slit to the beginning of the next (that is, the periodicity of the grating is d). We impose an incident plane wave of the form u inc (x, y, t) = cos (ωt + ky), where k = 2π/a and ω = k/c, where c is the wave speed. Periodic BCs at x = ±d/2 (determining the periodicity of the grating), and homogeneous Dirichlet BCs are imposed at the screen. We also test outflow boundary conditions in multiple dimensions, which are imposed at y = ±L y /2.   In Figure 10, we observe the time evolution of the incident plane wave passing through the aperture, and the resulting interference patterns as the diffracted wave propagates across the periodic boundaries. The outflow boundary conditions allow the waves to propagate outside the domain. While a rigorous analysis of the efficacy of our outflow BCs is the subject of future work, the results look quite reasonable, as no spurious reflections are seen at the artificial boundaries. 5.5. Point Sources and Outflow Boundary Conditions. As a final example, we illustrate some of the more interesting features of our solver. We launch two point sources (from points not on the mesh), and employ Dirichlet, periodic and outflow boundary conditions along the edges of the domain. As can be seen in Figure 11, the point sources propagate perfectly, despite not being placed at grid points.

Conclusion
In this paper we have presented a fast, A-stable, second order scheme for solving the wave equation. Using the method of lines transpose (MOL T ), the solution can be interpreted in the semi-discrete sense as a boundary integral solution, posed as a convolution against an exponential kernel. We have exploited this fact to develop a matrix-free, O(N ) fast spatial convolution algorithm, capable of embedding boundary points in a regular Cartesian mesh, without affecting the accuracy or stability.
In addition to demonstrating second order accuracy for wave propagation in a variety of non-Cartesian geometries, with time steps in excess of the usual CFL restriction, we have also developed a novel method for implementing outflow boundary conditions, as well as methods for launching waves, using soft sources, from points located at arbitrary locations (e.g., not located at a mesh point) inside the domain, which is of interest in particle-wave simulations such as those required in studying plasma. Several topics warrant further investigation. We are developing a domain decomposition approach to multi-core computing with our implicit wave solver based on the properties of the exponential kernel, where the subdomains require only pointwise data communication from adjoining edges. Future work will investigate the implementation of higher-order accurate boundary conditions for complex boundary geometries and outflow boundary conditions, and further, apply these methods to Maxwell's equations both in electromagnetic scattering and plasma physics problems.
Appendix A. Summary Table For Second-Order Wave Solver
We are predominantly interested in the case where S(x, t) consists of a large number of time dependent point sources. However, it is often the case that in electromagnetics problems, a soft source is prescribed to excite waves of a prescribed frequency, or range of frequencies, within the domain. A soft source is so named because, although incident fields are generated at a prescribed fixed spatial location, no scattered fields are generated.
The implementation of a soft source σ(t) at x = x s is accomplished by prescribing the source condition u(x s , t) = σ(t).
(48) However, it can be shown that if we set and insert it into the wave equation, then the soft source condition (48) is satisfied, and the solutions are equivalent. Thus, a soft source is nothing more than a point source, whose time-varying field is integrated by the wave equation. Upon convolving this source term with the Green's function according to (4), we find where the definition of α = β/(c∆t) has been utilized.
Remark 5. It is often the case that taking the analytical derivative σ (t n ) is to be avoided, for various reasons. In this case, any finite difference approximation which is of the desired order of accuracy can be substituted.
Likewise for general point sources, the corresponding form of the source term is Therefore, it suffices to consider delta functions both for the implementation of soft sources, as well as including time dependent point sources.