A unified integral equation scheme for doublyperiodic Laplace and Stokes boundary value problems in two dimensions
Abstract
We present a spectrallyaccurate scheme to turn a boundary integral formulation for an elliptic PDE on a single unit cell geometry into one for the fully periodic problem. Applications include computing the effective permeability of composite media (homogenization), and microfluidic chip design. Our basic idea is to exploit a small least squares solve to apply periodicity without ever handling periodic Green’s functions. We exhibit fast solvers for the twodimensional (2D) doublyperiodic Neumann Laplace problem (flow around insulators), and Stokes nonslip fluid flow problem, that for inclusions with smooth boundaries achieve 12digit accuracy, and can handle thousands of inclusions per unit cell. We split the infinite sum over the lattice of images into a directlysummed “near” part plus a small number of auxiliary sources which represent the (smooth) remaining “far” contribution. Applying physical boundary conditions on the unit cell walls gives an expanded linear system, which, after a rank1 or rank3 correction and a Schur complement, leaves a wellconditioned square system which can be solved iteratively using fast multipole acceleration plus a lowrank term. We are rather explicit about the consistency and nullspaces of both the continuous and discretized problems. The scheme is simple (no lattice sums, Ewald methods, nor particle meshes are required), allows adaptivity, and is essentially dimension and PDEindependent, so would generalize without fuss to 3D and to other nonoscillatory elliptic problems such as elastostatics. We incorporate recently developed spectral quadratures that accurately handle closetotouching geometries. We include many numerical examples, and provide a software implementation.
1 Introduction
Periodic boundary value problems (BVPs) arise frequently in engineering and the sciences, either when modeling the behavior of solid or fluid media with true periodic geometry, or when applying artificial periodic boundary conditions to simulate a representative domain of a random medium or particulate flow (often called a supercell or representative volume element simulation). The macroscopic response of a given microscopic periodic composite medium can often be summarized by an effective material property (e.g. a conductivity or permeability tensor), a fact placed on a rigorous footing by the field of homogenization (for a review see [14]). However, with the exception of layered media that vary in only one dimension, finding this tensor requires the numerical solution of “cell problems” [60, 64], namely BVPs in which the solution is periodic up to some additive constant expressing the macroscopic driving. Application areas span all of the major elliptic PDEs, including the Laplace equation (thermal/electrical conductivity, electrostatics and magnetostatics of composites [35, 38, 27, 12]); the Stokes equations (porous flow in periodic solids [18, 49, 73, 26], sedimentation [1], mobility [67], transport by cilia carpets [16], vesicle dynamics in microfluidic flows [57]); elastostatics (microstructured periodic or random composites [29, 36, 59, 62]); and the Helmholtz and Maxwell equations (phononic and photonic crystals, bandgap materials [41, 63]). In this work we focus on the first two (nonoscillatory) PDEs above, noting that the methods that we present also apply with minor changes to the oscillatory Helmholtz and Maxwell cases, at least up to moderate frequencies [5, 13, 55].
The accurate solution of periodic BVPs has remained a challenging problem at the forefront of analytical and computational progress for well over a century [9]. Modern simulations may demand large numbers of objects per unit cell, with arbitrary geometries, that at high volume fractions may approach arbitrarily close to each other [70]. In such regimes asymptotic methods based upon expansion in the inclusion size do not apply [18, 26]. Polydisperse suspensions or nonsmooth geometries may require spatially adaptive discretizations. In microfluidics, high aspect ratio and/or skew unit cells are needed, for instance in the optimal design of particle sorters. Furthermore, to obtain accurate average material properties for random media or suspensions via the Monte Carlo method, thousands of simulation runs may be needed [38]. A variety of numerical methods are used (see [60, Sec. 2.8]), including particular solutions (starting in 1892 with Rayleigh’s method for cylinders and spheres [68], and, more recently, in Stokes flow [69]), eigenfunction expansions [73], lattice Boltzmann and finite differencing [47], and finite element methods [29, 59]. However, for general geometries, it becomes very hard to achieve highorder accuracy with any of the above apart from finite elements, and the cost of meshing renders the latter unattractive when moving geometries are involved. Integral equation methods [48, 65, 44, 54] are natural, since the PDE has piecewiseconstant coefficients, and are very popular [49, 36, 27, 26, 62, 12, 1, 67]. By using potential theory to represent the solution field in terms of the convolution of a freespace Green’s function with an unknown “density” function living only on material boundaries, one vastly reduces the number of discretized unknowns. The linear system resulting by applying boundary conditions takes the form
\hb@xt@.01(1.1) 
where the matrix is the discretization of an integral operator. The latter is often of Fredholm secondkind, hence remains wellconditioned, independent of , so that iterative methods converge rapidly. Further advantages include their highorder or spectral accuracy (e.g. via Nyström quadratures), the ease of adaptivity on the boundary, and the existence of fast algorithms to apply with optimal cost, such as the fast multipole method (FMM) [28].
The traditional integral equation approach to periodic problems replaces the freespace Green’s function by a periodic Green’s function, so that the density is solved on only the geometry lying in a single unit cell [36, 27, 26, 62, 12]. There is an extensive literature on the evaluation of such periodic Green’s functions (e.g. in the Stokes case see [66, 72]); yet any such pointwise evaluation for each sourcetarget pair leads to complexity, which is unacceptable for large problems. There are two popular approaches to addressing this in a way compatible with fast algorithms:

Lattice sums. Noticing that the difference between the periodic and freespace Green’s function is a smooth PDE solution in the unit cell, one expands this in a particular solution basis (a cylindrical or spherical expansion); the resulting coefficients, which need to be computed only once for a given unit cell, are called lattice sums [7, 35, 27, 26]. They originate in the work of Rayleigh [68] and in the study of ionic crystals (both reviewed in [9, Ch. 2–3]). The FMM may then be periodized by combining the toplevel multipole expansion coefficients with the lattice sums to give a correction to the local expansion coefficients (in 2D this is a discrete convolution) which may be applied fast [28, Sec. 4.1] [62]. This method has been used for the doublyperiodic Laplace BVP by Greengard–Moura [27], and Stokes BVP by Greengard–Kropinski [26].

Particlemesh Ewald (PME) methods. These methods exploit Ewald’s realization [20] that, although both the spatial and the spectral (Fourier) sums for the periodic Green’s function in general converge slowly, there is an analytic splitting into spatial and spectral parts such that both terms converge superalgebraically. Numerically, the spatial part is now local, while the spectral part may be applied by “smearing” onto a uniform grid, using a pair of FFTs, and (as in the nonuniform FFT [19]) correcting for the smearing. The result is a fast algorithm; for a review see [15]. Recently this has been improved by Lindbo–Tornberg to achieve overall spectral accuracy in the Laplace [52] and Stokes [50] settings, with applications to 3D fluid suspensions [1].
The scheme that we present is based purely on freespace Green’s functions, and has distinct advantages over the above two periodization approaches. 1) Only physicallymeaningful boundary conditions and compatibility conditions are used. Aside from conceptual and algebraic simplification, this also removes the need for (often mysterious and subtle [9]) choices of the values of various conditionally or even nonconvergent lattice sums based on physical arguments [18, 28, 36, 27] (e.g., six such choices are needed in [26]). When it comes to skewed unit cells combined with applied pressure driving, such handpicked choices in lattice sum methods become far from obvious. 2) Freespace FMM codes may be inserted without modification, and fast direct solvers (based on hierarchical compression of matrix inverses) [31, 23] require only slight modification [22, 57]. 3) When multiscale spatial features are present in the problem—commonly observed in practical applications such as polydisperse suspensions and complex microfludic chip geometries—PME codes become inefficient because of their need for uniform grids. In contrast, our algorithm retains the fast nature of the adaptive FMM in such cases. 4) In contrast to lattice sums and Ewald methods, which intrinsically rely on radiallysymmetric expansions and kernels, our scheme can handle reasonably high aspect ratio unit cells with little penalty. 5) Since our Stokes formulation does not rely on complex variables (e.g. Sherman–Lauricella [26]), the whole scheme generalizes to 3D with no conceptual changes [56]. A disadvantage of our scheme is that the prefactor may be slightly worse than Ewald methods, due to the direct summation of neighboring images.
The basic idea is simple, and applies to integral representations for a variety of elliptic PDEs. Let be a (generally trapezoidal) unit cell “box”, and let be the material boundary which is to be periodized ( may even intersect the walls of ). Applying periodic boundary conditions is equivalent to summing the freespace representation on over the infinite lattice, as in Fig. LABEL:f:geom(a). Our scheme sums this potential representation over only the nearest neighbor images of (see Fig. LABEL:f:geom(b)), but adds a small auxiliary basis of smooth PDE solutions in , with coefficient vector , that efficiently represents the (distant) contribution of the rest of the infinite image lattice. For this auxiliary basis we use point sources (Fig. LABEL:f:geom(c)); around are needed, a small number which, crucially, is independent of and the boundary complexity. We apply the usual (homogeneous) boundary condition on , which forms the first block row of a linear system. We impose the desired physical periodicity as auxiliary conditions on the Cauchy data differences between wall pairs – and – (see Fig. LABEL:f:geom(a)), giving a second block row. The result is a block “extended linear system” (ELS),
\hb@xt@.01(1.2) 
Here accounts for the applied macroscopic thermal or pressure gradient in the form of prescribed jumps across one unit cell. Fig. LABEL:f:geom(b)–(c) sketches the interactions described by the four operators , , , and . This linear system is generally rectangular and highly illconditioned, but when solved in a backwardstable leastsquares sense can result in accuracies close to machine precision.
Three main routes to a solution of (LABEL:ELS) are clear. (a) The simplest is to use a dense direct leastsquares solve (e.g., via QR), but the cost becomes impractical for large problems with . Instead, to create fast algorithms one exploits the fact that the numbers of auxiliary rows and columns in (LABEL:ELS) are both small, as follows. (b) One may attempt to eliminate to get the Schur complement square linear system involving , the periodized version of , of the form
\hb@xt@.01(1.3) 
where is a pseudoinverse of (see Section LABEL:s:schur). As we will see, (LABEL:Aper) can be wellconditioned when (LABEL:A) is, and can be solved iteratively, by using the FMM to apply while applying the second term in its lowrank factored form . (c) One may instead eliminate by forming, then applying, a compressed representation of via a fast direct solver; this has proven useful for the case of fixed periodic geometry with multiple boundary data vectors [22, 57]. Both methods (b) and (c) can achieve complexity.
This paper explores route (b). A key contribution is overcoming the stumbling block that, for many standard PDEs and integral representations, as written in (LABEL:Aper) does not in fact exist. Intuitively, this is due to divergence in the sum of the Green’s function over the lattice. For instance, the sum of clearly diverges, and thus the singlelayer representation for Laplace’s equation, which we use in Section LABEL:s:lap, cannot be periodized unless the integral of density (total charge) vanishes. This manifests itself in the 2nd block row of (LABEL:ELS): the range of contains vectors that are not in the range of , thus numerically blows up. After studying in Sections LABEL:s:lapempty and LABEL:s:stoempty the consistency conditions for the linear system involving (the “empty unit cell” problem), we propose rank1 (for Neumann Laplace) and rank3 (for noslip Stokes) corrections to the ELS that allow the Schur complement to be taken, and, moreover that remove the nullspace associated with the physical BVP. We justify these schemes rigorously, and test them numerically.
Although the idea of imposing periodicity via extra linear conditions has aided certain electromagnetics solvers for some decades [32, 75], we believe that this idea was first combined with integral equations by the first author and Greengard in [5], where controlled a Helmholtz local expansion. Since then it has become popular for a variety of PDE [5, 22, 30, 13, 12, 57, 55]. From [5] we also inherit the split into near and far images, and cancellations in that allow a rapidly convergent scheme even when intersects unit cell walls. The use of point sources as a particular solution basis that is efficient for smooth solutions is known as the “method of fundamental solutions” (MFS) [8, 21], “method of auxiliary sources” [45], “charge simulation method” [43], or, in fast solvers, “equivalent source” [11] or “proxy” [58] representations. This is also used in the recent 3D periodization scheme of Gumerov–Duraiswami [30]. Finally, the lowrank perturbations which enlarge the range of are inspired by the lowrank perturbation methods for singular square systems of Sifuentes et al. [71].
Here is a guide to the rest of this paper. In Section LABEL:s:lap we present the periodic Neumann Laplace BVP, possibly the simplest application of the scheme. This first requires understanding and numerically solving the “empty unit cell” subproblem, which we do in Section LABEL:s:lapempty. The integral formulation, discretization, and direct numerical solution of the ELS (LABEL:ELS) follows in Section LABEL:s:lapbie. A general scheme to stably take the Schur complement is given in Section LABEL:s:schur, along with a scheme to remove the physical nullspace specific to our representation. The latter is tested numerically. Section LABEL:s:multi shows how the FMM and closeevaluation quadratures are incorporated, and applies them to large problems with thousands of inclusions. Section LABEL:s:kappa defines the effective conductivity tensor and shows how to evaluate it efficiently using a pair of BVP solutions. In Section LABEL:s:sto we move to periodic Dirichlet (noslip) Stokes flow, and show how its periodizing scheme closely parallels the Laplace version. In fact, the only differences are new consistency conditions for the empty subproblem (Section LABEL:s:stoempty), the use of a combinedfield formulation (Section LABEL:s:stobie), and the need for a rankthree perturbation for a stable Schur complement (Section LABEL:s:stosch). Experiments on the drag of a square array of discs, and on largescale flow problems, are performed in Section LABEL:s:stonum. We discuss generalizations and conclude in Section LABEL:s:conc.
Our aim is to illustrate, with two example BVPs, a unified route to a wellconditioned periodization compatible with fast algorithms. We take care to state the consistency conditions and nullspaces of the full and empty problems, this being the main area of problemspecific variation. We believe that the approach adapts simply to other boundary conditions and PDEs, once the corresponding consistency conditions and nullspaces are laid out. Although general aspect ratios and skew unit cells are amongst our motivations, for clarity we stick to the unit square; the generalization is straightforward.
Remark 1.1 (Software)
We maintain documented MATLAB codes implementing the basic methods in this paper at http://github.com/ahbarnett/BIE2D
2 The Neumann Laplace case
We now present the heat/electrical conduction problem in the exterior of a periodic lattice of insulating inclusions (corresponding to in [27]). For simplicity we first assume a single inclusion per unit cell. Let and be vectors defining the lattice in ; we work with the unit square so that and . Let represent the infinite lattice of inclusions. The scalar , representing electric potential or temperature, solves the BVP
\hb@xt@.01(2.1)  
\hb@xt@.01(2.2)  
\hb@xt@.01(2.3)  
\hb@xt@.01(2.4) 
ie, is harmonic, has zero flux on inclusion boundaries, and is periodic up to a given pair of constants which encode the external (macroscopic) driving. We use the abbreviation , where is the unit normal on the relevant evaluation curve.
Proposition 2.1
For each the solution to (LABEL:lap)–(LABEL:p2) is unique up to an additive constant.
Proof. As usual, one considers the homogeneous BVP arising when is the difference of two solutions. Let be any unit cell (tiling domain) containing , then using Green’s first identity,
The first boundary term cancels by periodicity, the second by (LABEL:bc), hence .
To solve the periodic BVP we first reexpress it as a BVP on a single unit cell with coupled boundary values on the four walls comprising its boundary ; see Fig. LABEL:f:geom(a). For simplicity we assume for now that the square can be chosen such that does not intersect any of the walls; this restriction will later be lifted (see Remark LABEL:r:intersect). We use the notation to mean the restriction of to the wall , and for its normal derivative using the normal on (note that this points to the right, as shown in the figure.) Consider the following reformulated BVP:
\hb@xt@.01(2.5)  
\hb@xt@.01(2.6)  
\hb@xt@.01(2.7)  
\hb@xt@.01(2.8)  
\hb@xt@.01(2.9)  
\hb@xt@.01(2.10) 
Clearly any solution to (LABEL:lap)–(LABEL:p2) also satisfies (LABEL:lap2)–(LABEL:j2). Because of the unique continuation of Cauchy data as a solution to the 2ndorder PDE, the converse holds, thus the two BVPs are equivalent. We define the discrepancy [5] of a solution as the stack of the four functions on the lefthand side of (LABEL:u1)–(LABEL:j2).
2.1 The empty unit cell discrepancy BVP and its numerical solution
We first analyze, then solve numerically, an important subproblem which we call the “empty unit cell BVP.” We seek a harmonic function matching a given discrepancy (i.e., a stack of four functions defined on the walls , , , , respectively). That is,
\hb@xt@.01(2.11)  
\hb@xt@.01(2.12)  
\hb@xt@.01(2.13)  
\hb@xt@.01(2.14)  
\hb@xt@.01(2.15) 
We now give the consistency condition and nullspace for this BVP, which shows that it behaves essentially like a square linear system with nullity 1.
Proposition 2.2
A solution to (LABEL:lapv)–(LABEL:g4) exists if and only if , and is then unique up to a constant.
Proof. The zeroflux condition holds for harmonic functions. Writing as the union of the four walls, with their normal senses, gives the sum of the integrals of (LABEL:g2) and (LABEL:g4). Uniqueness follows from the method of proof of Prop. LABEL:p:lapnull.
We now describe a numerical solution method for this BVP that is accurate for a certain class of data , namely those for which the solution may be continued as a regular solution to Laplace’s equation into a large neighborhood of , and hence each of is analytic. It will turn out that this class is sufficient for our periodic scheme (essentially because will only have to represent distant image contributions, as shown in Fig. LABEL:f:geom(c)). Let be centered at the origin. Recalling the fundamental solution to Laplace’s equation,
\hb@xt@.01(2.16) 
we approximate the solution in by a linear combination of such solutions with source points lying uniformly on a circle of radius , where is the maximum radius of the unit cell. That is, for ,
\hb@xt@.01(2.17) 
with unknown coefficient vector . As discussed in the introduction, this is known as the MFS. Each basis function is a particular solution to (LABEL:lapv) in , hence only the boundary conditions need enforcing (as in Rayleigh’s original method [68]).
This MFS representation is complete for harmonic functions. More precisely, for in a suitable class, it is capable of exponential accuracy uniformly in . In particular, since is contained within the ball , we may apply known convergence results to get the following.
Theorem 2.3
Let extend as a regular harmonic function throughout the closed ball of radius . Let the fixed proxy radius satisfy . For each let be a proxy basis set as in (LABEL:mfs). Then there is a sequence of coefficient vectors , one vector for each , and a constant dependent only on , such that
In addition, the vectors may be chosen so that the sequence , , is bounded.
This exponential convergence—with rate controlled by the distance to the nearest singularity in —was derived by Katsurada [42, Rmk. 2.2] in the context of collocation (point matching) for a Dirichlet BVP on the boundary of the disc of radius , which is enough to guarantee the existence of such a coefficient vector sequence. The boundedness of the coefficient norms has been known in the context of Helmholtz scattering in the Russian literature for some time (see [46, Sec. 3.4] and references within, and [17, Thm. 2.4]). We do not know of a reference stating this for the Laplace case, but note that the proof is identical to that in [4, Thm. 6], by using equation (13) from that paper. The restriction is crucial for high accuracy, since when , although convergence occurs, blows up, causing accuracy loss due to catastrophic cancellation [4, Thm. 7].
Remark 2.4
Intuitively, the restriction arises because the proxy basis may be viewed as a periodic trapezoid rule quadrature approximation to a singlelayer potential on the circle , a representation which is incomplete when : it cannot represent the constant function [74, Rmk. 1].
To enforce boundary conditions, let , , , be two sets of collocation points, on the left and bottom wall respectively; we use Gauss–Legendre nodes. Enforcing (LABEL:g1) between collocation points on the left and right walls then gives
\hb@xt@.01(2.18) 
Continuing in this way, the full set of discrepancy conditions (LABEL:g1)–(LABEL:g4) gives the linear system
\hb@xt@.01(2.19) 
where stacks the four vectors of discrepancies sampled at the collocation points, while the matrix consists of four block rows, each of size . By reading from (LABEL:Q1), one sees that has elements . Analogously, , , and finally .
We solve the (generally rectangular) system (LABEL:Qsys) in the least squares sense, which corresponds to minimizing a weighted norm of the discrepancy error (we find in practice that there is no advantage to incorporating the squareroot of the weights associated with the collocation nodes). As is well known in the MFS community [42, 21], becomes exponentially illconditioned as grows; intuitively this follows from the exponential decay of the singular values of the singlelayer operator from the proxy radius to the circle radius containing . Thus, a least squares solver is needed which can handle rankdeficiency stably, i.e. return a smallnorm solution when one exists. In this case the illconditioning causes no loss of accuracy, and, even though there is instability in the vector , for the evaluation of errors close to machine precision () are achieved [4].
Remark 2.5
For this and subsequent direct solves we use MATLAB’s linsolve (with the option RECT=true to prevent the less accurate LU from being used in the square case), which uses columnpivoted QR to find the socalled “basic” solution [61, Sec. 5.7] having at most nonzero entries, where is the numerical rank; this is close to having minimum norm [25, Sec. 5.5].
We illustrate this with a simple numerical test, in which is the discrepancy of a known harmonic function of typical magnitude and with sufficiently distant singularity. Once is filled and (LABEL:Qsys) solved, the numerical solution is evaluated via (LABEL:mfs), and the maximum error at 100 random target points in is taken (after removing an overall constant error, expected from Prop. LABEL:p:lapemptynull). Fig. LABEL:f:Qconv(a) shows exponential convergence in this error vs the number of wall nodes . Fig. LABEL:f:Qconv(b) shows exponential convergence with respect to , the number of proxy points, with a rate slightly exceeding the predicted rate. It is clear that whenever and the norm remains , and around 15 digits of accuracy result.
The decaying lowest few singular values of are also shown in panel (b): in particular there is one singular value decaying faster than all others. It is easy to verify that this corresponds to the nullspace of the BVP (Prop. LABEL:p:lapemptynull); indeed we have tested that its right singular vector is approximately constant, and generates via (LABEL:mfs) the constant function in to within relative error. Likewise, the consistency condition in Prop. LABEL:p:lapemptynull manifests itself in : let be the vector that applies the discretization of this consistency condition to a vector , namely
\hb@xt@.01(2.20) 
where semicolons indicate vertical stacking, are the vectors of weights corresponding to the collocation nodes on , , and is the zero vector. Then we expect that
\hb@xt@.01(2.21) 
and indeed observe numerically that once . In summary, although the matrix is generally rectangular and illconditioned, it also inherits both aspects of the unit nullity of the empty BVP that it discretizes.
Remark 2.6
A different scheme is possible in which would be square, and (modulo a nullity of one, as above) wellconditioned, based on a “tictactoe” set of layer potentials (see [5, Sec. 4.2] in the Helmholtz case). However, we recommend the above proxy point version, since i) the matrix is so small that handling its illconditioning is very cheap, ii) the tictactoe scheme demands closeevaluation quadratures for its layer potentials, and iii) the tictactoe scheme is more complicated.
2.2 Extended linear system for the conduction problem
We now treat the above empty BVP solution scheme as a component in a scheme for the periodic BVP (LABEL:lap2)–(LABEL:j2). Simply put, we take standard potential theory for the Laplace equation [44, Ch. 6], and augment this by enforcing periodic boundary conditions. Given a density function on the inclusion boundary , recall that the singlelayer potential, evaluated at a target point , is defined by
\hb@xt@.01(2.22) 
Using to indicate the outward normal at , this potential obeys the jump relation [44, Thm. 6.18]
\hb@xt@.01(2.23) 
where denotes the usual transposed doublelayer operator from a general source curve to a target curve , defined by
\hb@xt@.01(2.24) 
The integral implied by the selfinteraction operator is to be interpreted in the principal value sense.
Our representation for the solution sums the singlelayer potential over the copies closest to the origin, and adds an auxiliary basis as in (LABEL:mfs),
\hb@xt@.01(2.25) 
Our unknowns are the density and auxiliary vector . Substituting (LABEL:urep) into the Neumann boundary condition (LABEL:bc2), and using the exterior jump relation on the central copy () only, gives our first block row,
\hb@xt@.01(2.26) 
where, as before, the “near” superscript denotes summation over source images as in (LABEL:urep).
The second block row arises as follows. Consider the substitution of (LABEL:urep) into the first discrepancy equation (LABEL:u1): there are nine source copies, each of which interacts with the and walls, giving 18 terms. However, the effect of the rightmost six sources on is cancelled by the effect of the leftmost six sources on , leaving only six terms, as in [5, Fig. 4(a)–(b)]. All of these surviving terms involve distant interactions (the distances exceed one period if is contained in ). Similar cancellations occur in the remaining three equations (LABEL:j1)–(LABEL:j2). The resulting four subblocks are
\hb@xt@.01(2.27)  
\hb@xt@.01(2.28)  
\hb@xt@.01(2.29)  
\hb@xt@.01(2.30) 
Remark 2.7 (Wall intersection)
The cancellation of all near interactions in (LABEL:row2a)–(LABEL:row2d) is due to the neighbor summation in the representation (LABEL:urep). Furthermore this cancellation allows an accurate solution even when intersects , without the need for specialized quadratures, as long as remains far from the boundary of the enclosing unit cell block. Informally, the unit cell walls are “invisible” to the inclusions. With an elongated for which the last condition does not hold, with a little bookkeeping one could split into pieces that, after lattice translations, lie far from the boundary of the enclosing unit cell block; we leave this last case for future work.
(LABEL:bie)–(LABEL:row2d) form a set of coupled integralalgebraic equations, where the only discrete aspect is that of the proxy points.^{1}^{1}1It would also be possible to write a purely continuous version by replacing the proxy circle by a singlelayer potential; however, sometimes an intrinsically discrete basis is useful [5, 13]. It is natural to ask how the unit nullity of the BVP manifests itself in the solution space for the pair . Are there pairs with no effect on , enlarging the nullspace? It turns out that the answer is no, and that the onedimensional nullspace (constant functions) is spanned purely by , as a little potential theory now shows.
Lemma 2.8
In the solution to (LABEL:bie)–(LABEL:row2d), is unique.
Proof. Let be the difference between any two solutions to (LABEL:bie)–(LABEL:row2d). Let be the representation (LABEL:urep) using this , both in but also inside . Then by construction is a solution to the homogeneous BVP (LABEL:lap2)–(LABEL:j2), i.e. with , thus by Prop. LABEL:p:lapnull, is constant in . Thus by the continuity of the singlelayer potential [44, Theorem 6.18], the interior limit of on is constant. However, is harmonic in , so is constant in . By (LABEL:JR), , but we have just shown that both and vanish, so .
2.2.1 Discretization of the ELS
We discretize (LABEL:bie), the first row of the system, using a set of quadrature nodes on and weights such that
holds to high accuracy for smooth functions . In practice, when is parametrized by a periodic function , , then using the periodic trapezoid rule in gives and . Since has a smooth kernel, we apply Nyström discretization [44, Sec. 12.2] to the integral equation (LABEL:bie) using these nodes, to get our first block row
\hb@xt@.01(2.31) 
where is a vector of density values at the nodes, and has entries
\hb@xt@.01(2.32) 
where is the Kronecker delta. Here, for entries the standard diagonal limit of the kernel is needed, where is the signed curvature at . The matrix has entries .
For the second block row (LABEL:row2a)–(LABEL:row2d) we use the above quadrature for the source locations of the operators, and enforce the four equations on the wall collocation nodes , , , , , to get
\hb@xt@.01(2.33) 
with the macroscopic driving encoded by the righthand side vector
\hb@xt@.01(2.34) 
where is the vector of ones. The matrix is precisely as in (LABEL:Qsys). Rather than list formulae for all four blocks in , we have , with the others filled analogously. Stacking the two block rows (LABEL:row1) and (LABEL:row2) gives the by extended linear system (LABEL:ELS), which we emphasize is just a standard discretization of the BVP conditions (LABEL:bc2)–(LABEL:j2).
For small problems ( less than a few thousand), the ELS is most simply solved by standard dense direct methods; for larger problems a Schur complement must be taken in order to solve iteratively, as presented in Section LABEL:s:schur. First we perform a numerical test of the direct method.
2.3 Numerical tests using direct solution of the ELS
Example 1. We define a smooth “worm” shaped inclusion which crosses the unit cell walls and by . The solution of the periodic Neumann Laplace BVP with external driving is shown in Fig. LABEL:f:ELS(a). Here is evaluated via (LABEL:urep) both inside (where it is accurate), and, to show the nature of the representation, out to the proxy circle (where it is certainly inaccurate). Recall that the proxy sources must represent the layer potentials due to the infinite lattice of copies which excludes the central block. Note that two tips of copies in this set slightly penetrate the proxy circle, violating the condition in Thm. LABEL:t:Qconv that the function the proxy sources represent be analytic in the closed ball of radius . We find in practice that such slight geometric violations do not lead to problematic growth in , but that larger violations can induce a large which limits achievable accuracy. Panel (c) shows the convergence of errors in (at the pair of points shown) to their numericallyconverged value, fixing converged values for and . Convergence to around 13 digits (for , ) is apparent, and the solution time is 0.03 seconds.
As an independent verification of the method, we construct a known solution to a slightly generalized version of (LABEL:lap2)–(LABEL:j2), where (LABEL:bc2) is replaced by inhomogeneous data and a general discrepancy is allowed. We choose the known solution
where the central dipole has direction , and location chosen inside and far from its boundary (so that the induced data is smooth). The grid size of (some of which is shown by symbols in Fig. LABEL:f:ELS(a)) is chosen so that is sufficiently smooth (which requires at least ), and so that the periodizing part is nontrivial. From the righthand side functions and are then evaluated at nodes, the ELS solved directly, the numerical solution (LABEL:urep) evaluated, and the difference at two points compared to its known exact value. The resulting convergence to 13 digits is shown in Fig. LABEL:f:ELS(c). The convergence rate is slower than before, due to the unavoidable closeness of the dipole source to .
Finally, the grey lines in panel (e) show that, as with , there is one singular value of which is much smaller than the others, reflecting the unit nullity of the underlying BVP (Prop. LABEL:p:lapnull). Also apparent is the fact that, despite the exponential illconditioning, the solution norm remains bounded once convergence has occurred.
2.4 Schur complement system and its iterative solution
When is large, solving the illconditioned rectangular ELS (LABEL:ELS) is impractical. We would like to use a Schur complement in the style of (LABEL:Aper) to create an equivalent system, which we do in Sec. LABEL:s:equiv. Furthermore, in order to use Krylov subspace iterative methods with known convergence rates, we would like to remove the nullspace to make this wellconditioned, which we do in Sec. LABEL:s:wellcond. We will do both these tasks via lowrank perturbation of certain blocks of (LABEL:Aper) before applying the Schur complement.
In what follows, is the pseudoinverse [25, Sec. 5.5] of , i.e. the linear map that recovers a smallnorm solution (if one exists) to via . The obstacle to using (LABEL:Aper) as written is that inherits a consistency condition from the empty BVP so that has one smooth vector in its left nullspace (; see (LABEL:wQ)). However the range of does not respect this condition, thus has a huge 2norm (which we find numerically is at least ). We first need to show that the range of a rankdeficient matrix may be enlarged by a rank perturbation, a rectangular version of results about singular square matrices in [71].
Lemma 2.9
Let have a dimensional nullspace. Let have fullrank projection onto , i.e. if has columns forming a basis for then is invertible. Let be arbitrary. Then , i.e. the range now includes that of .
Proof. We need to check that has a solution for all given pairs , . Recalling that is invertible, by substitution one may check that is an explicit such solution.
Remark 2.10
With additional conditions , and that has fullrank projection onto a part of (i.e. if has columns forming a basis for a dimensional subspace of , then is invertible), we get incidentally that has trivial nullspace (generalizing [71, Sec. 3]). The proof is as follows. Let solve the homogeneous equation . Then , but is invertible, so . Thus the homogeneous equation becomes , which means there is an such that . Thus , but is invertible, so that , so .
We also need the fact that a blockcolumn operation allows to be perturbed as above while changing the ELS solution space in a known way. The proof is simple to check.
Proposition 2.11
Let , , and be matrices, and let be a matrix with as many rows as has columns and as many columns as has rows. Then the pair solves the block system
\hb@xt@.01(2.35) 
if and only if the pair , with , solves the original ELS (LABEL:ELS).
2.4.1 An equivalent square system preserving the nullspace
Armed with the above, a method to “fold” the Neumann Laplace ELS into an equivalent square system is as follows:

Set the proxy coefficient vector , and the discrete density vector . Create lowrank perturbed matrix blocks and .

Solve for the vector in the Schur complement linear system
\hb@xt@.01(2.36) More precisely, since (due to the numerical illconditioning in ) multiplying by would lose accuracy due to rounding error, instead solve the small systems for , and for , heeding Remark LABEL:r:linsolve. From them, build and , which are respectively the system matrix and righthand side for (LABEL:Gschur). This large square system may then be solved iteratively (see section LABEL:s:multi).

Recover the proxy coefficients via .

Recover the density via .
Note that in Step 1 the prefactor leads to correct quadrature scaling, so that has similar 2norm to the other matrices (also recommended in [71]).
Theorem 2.12
Let encode the driving as in (LABEL:gg). Then for all sufficiently large , , and , any pair produced by the above procedure performed in exact arithmetic solves the original ELS (LABEL:ELS).
Proof. By rotational invariance, a constant singlelayer density on a circle generates a constant potential inside, and inserting into (LABEL:mfs) gives the periodic trapezoid quadrature approximation to such a potential, thus generates discrepancy near zero (in fact exponentially small in ). Thus for all sufficiently large , is not orthogonal to . Applying Lemma LABEL:l:onesmatrix with gives that the range of includes the discrepancy vector produced by the constant density . Thus, to show that the range of includes all smooth vectors, i.e. that it does not have the consistency condition (LABEL:wQ), one needs to check that which is done in Lemma LABEL:l:onesdens below. Thus the systems and are consistent for any and , so that the Schur complement is well defined. Finally, Prop. LABEL:p:gary, using the rank1 matrix , insures that step 4 recovers a solution to (LABEL:ELS).
For the missing technical lemma, we first need a form of Gauss’ law, stating that for any density on a curve the singlelayer potential generates flux equal to the total charge, i.e.
\hb@xt@.01(2.37) 
where is the boundary of some open domain containing . This may be proved by combining the jump relation (from (LABEL:JR)) with the fact that, since is a Laplace solution in , taken over the boundaries of and of .
Lemma 2.13
Let the matrix be defined as in Section LABEL:s:lapbie, and be as in (LABEL:w). Then, for all and sufficiently large, .
Proof. Let be the discrepancy of the potential for density . As used in the proof of Prop. LABEL:p:lapemptynull, the flux (left hand side of (LABEL:Gauss)) out of the unit cell equals , which by (LABEL:Gauss) (and noting that only one of the nine source terms in lies within ) equals the perimeter . What we seek is the discretization of this statement about the PDE. If is the discrepancy of sampled at the wall nodes, then its quadrature approximation is , and so (as discussed above (LABEL:w)) the quadrature approximation to is . Thus, as and tend to infinity, the latter converges to .
Theorem LABEL:t:gary justifies rigorously one procedure to create a square system equivalent to the ELS (LABEL:ELS). However, by the equivalence in Prop. LABEL:p:gary, the system matrix at the heart of the procedure is singular, as it inherits the unit nullity of the ELS, which itself derives from the unit nullity of the Laplace Neumann BVP. Since the convergence of iterative solvers for singular matrices is a subtle matter [10], this motivates the following improved variant which removes the nullspace.
2.4.2 A wellconditioned square system
The following simpler variant creates a nonsingular square system from the ELS; its proof is more subtle. It is what we recommend for the iterative solution of the periodic Neumann Laplace problem, and test numerically:

Set the constant proxy coefficient vector , the constant discrete density vector , and .

Solve for the vector in the Schur complement linear system
\hb@xt@.01(2.38) where, as before, one solves the small systems and , to build the large system matrix and righthand side for (LABEL:pschur), which may then be solved iteratively (see section LABEL:s:multi).

Recover the proxy coefficients via