An integral equation based numerical method for the forced heat equation on complex domains

Integral equation based numerical methods are directly applicable to homogeneous elliptic PDEs, and offer the ability to solve these with high accuracy and speed on complex domains. In this paper, extensions to problems with inhomogeneous source terms and time dependent PDEs, such as the heat equation, have been introduced. One such approach for the heat equation is to first discretise in time, and in each time-step solve a so-called modified Helmholtz equation with a parameter depending on the time step size. The modified Helmholtz equation is then split into two parts: a homogeneous part solved with a boundary integral method and a particular part, where the solution is obtained by evaluating a volume potential over the inhomogeneous source term over a simple domain. In this work, we introduce two components which are critical for the success of this approach: a method to efficiently compute a high-regularity extension of a function outside the domain where it is defined, and a special quadrature method to accurately evaluate singular and nearly singular integrals in the integral formulation of the modified Helmholtz equation for all time step sizes.


Introduction
In this paper we present a highly accurate numerical method for solving the forced isotropic heat equation with Dirichlet data on complex multiple connected domains in two dimensions. We adapt the solution methodology introduced by Kropinski and Quaife in [1], but extend and generalise their work to allow for solution of a wider class of problems with improved discretisation in time and uniform accuracy all the way up to the boundary. First, the heat equation is discretized in time with an implicit treatment of the diffusion term, an approach that is sometimes referred to as Rothe's method [2,3] or elliptic marching. This results in a sequence of modified Helmholtz equations, also known as the linearised Poisson-Boltzmann equation, to be solved at each time step. Doing so advances the solution to the parabolic heat equation in time. A relaxed definition of the modified Helmholtz equation reads α 2 u − ∆u = f , with α 2 inversely proportional to the time step. Utilising the linearity, this equation is further split into two: one that finds a particular solution for the specific right hand side without enforcing the boundary conditions, and a homogeneous problem that ensures that the sum of the two solutions solves the original problem. The homogeneous problem is solved with a boundary integral method with a panel-based Nyström quadrature scheme, as introduced in [4] by Kropinski and Quaife. The particular solution is written as a volume potential with the free space Green's function for the modified Helmholtz equation, also known as the Yukawa-or screened Columb potential. To avoid constructing quadrature methods for the evaluation of this volume potential over complex domains, an extension of the right hand side f is introduced, allowing for integration over a simple rectangular domain.
In [1] the authors Kropinski and Quaife demonstrated the potential of developing an efficient and accurate general boundary integral solver for the heat equation on complex domains. Moreover, they list the major remaining issues that require further investigation. At that time only examples for which a continuous extension of f could be constructed by hand was considered, thus excluding complex geometries and general data. Another impediment was the loss of accuracy for evaluating layer potentials close to their sources. Their solution was to over-resolve the boundary, but the loss of accuracy is still significant as a target point approaches the boundary. In this paper we introduce the following developments: • High order adaptive methods for time evolution.
• A method to efficiently compute a high-regularity extension of a function f to an enclosing and geometrically simple domain, given only its values at discrete locations in Ω.
• A special purpose quadrature method to avoid loss of accuracy when evaluating layer potentials close to the boundary and the kernel becomes nearly singular.
Two main groups of semi-implicit time stepping methods are Runge-Kutta methods [5] and spectral deferred correction methods [6,7,8]. We use the former to obtain an adaptive scheme, but the approach we propose is general with respect to the choice of semi-implicit time stepper.
It is not a simple problem to construct a high regularity extension of a function, for which only its values are known in discrete points inside the original domain Ω. In [1], Kropinski and Quaife considered only examples for which a continuous extension could be constructed by hand. We use a partition of unity extension technique (PUX) by Fryklund, Lehto, and Tornberg in [9]. They solve the Poisson equation with the above-mentioned split into a particular and an homogeneous problem. We now use this method for function extension in the context of the modified Helmholtz equation with excellent results and can hence increase the class of solvable problems as compared to [1]. An alternative approach for function extension is given in [10], where the function to be extended outside of Ω sets the boundary Dirichlet data on ∂Ω for the Laplace equation in R 2 \ Ω. The solution to this problem is computed with an integral equation based method, and defines a continuous function extension. See [9] and the references therein for other extension techniques, such as Fourier continuation methods or extending the unknown solution or solution from previous time step [11,12,13].
When evaluating a layer potential close to a boundary, the kernel becomes nearly singular. A well known challenge with boundary integral based methods is accurate numerical integration of singular (for evaluation on the boundary) and nearly singular kernels. The comparative study [14] complemented with [15] give an overview of state of the art methods. The latter includes panel-based explicit kernel-split schemes with product integration, pioneered by Helsing and Ojala [16] for the Laplace equation. This methodology is applicable to a large class of linear elliptic PDEs, and achieves excellent results also for e.g. the Helmholtz [15] and Stokes equations [17]. However, for the modified Helmholtz equation product integration may fail altogether for sufficiently large α, i.e. for small time steps in our setting. The quadrature rule will in this case require an unfeasibly high resolution of the boundary, which is not motivated by the geometry nor the resolution requirement for the layer density. This spurred the development of a quadrature scheme to solve this problem. The Yukawa potential decays as exp (−α), and the kernel becomes more localised as α increases. In this process, the product integration requires an increasing amount of upsampling, but only over a decreasing interval, and hence only local upsampling is needed. In a separate paper [18], we present an adaptive quadrature scheme in the spirit of [19] that lifts the previous restriction on α.
A parallel development of a boundary integral based solver for the heat equation is based on direct approximation of the heat kernel, thus avoiding discretisation of the differential operator with respect to time. In the initial work [20] it was observed that to achieve the desired accuracy for domains with high curvature the time step must be considerably smaller than the formal rate of convergence would suggest. The authors refer to this as geometrically induced stiffness. In recent work towards solving the heat equation with said method Wang et al. has developed a hybrid method that allows for evaluation of the boundary and volume potentials including the space-time heat kernel without the constraints from geometric stiffness [21].
Efforts to solve the heat equation with boundary integral equation based techniques are not only motivated by that specific task. Surely, there are other methods to solve the heat equation on a complex domain, such as finite element methods. However, the algorithmic development in these efforts is essential to increase the applicability of The boundaries are denoted Γ n , n = 0, . . . , N Γ . The outer boundary is Γ 0 and the outward directed normal is denoted by ν.
integral equation based numerical methods which sport several attractive features, including that complex geometry naturally enters the problem and generation of an unstructured mesh is redundant, ill-conditioning associated with discretising the operators is avoided, high accuracy can be attained, and boundary data and far field conditions are simple to incorporate. Developments for the heat equation are also related to extension from Stokes to Navier-Stokes equations.
The focus of this paper is on the heat equation. However, fast integral equations for the modified Helmholtz equation are of interest for the many applications that equation applies to. These include, but are not limited to: electrostatic interactions in protein and related biological functions, macroscopic electrostatics, DebyeHuckel theory for dilute electrolytes, water wave problems, in the linearisation of the PoissonBoltzmann equation and approximation of surfaces [22,23,24,25,26,27]. Consequently, there is active research on solution methods and analysis thereof for the modified Helmholtz equation. In [28] the method of fundamental solution is used, while in [29] it is solved by plane wave functions.

Overview of the paper
The mathematical problem is formulated in Section 2, both for the heat equation and the modified Helmholtz equation. Section 3 contains the numerical methods for solving the homogeneous problem and the particular problem for the modified Helmholtz equation, including an introduction to PUX. It is assumed that the heat equation as been appropriately discretised in time. Thereafter we present the numerical results in Section 4, for the modified Helmholtz equation, the heat equation and a reaction-diffusion type problem. Finally we present our conclusions and an outlook in Section 5. See Appendix A.1 for instructions on how IMEX Runge-Kutta methods reduce the heat equation to a sequence of modified Helmholtz equations. There are simple examples, Butcher tableaus and a note on adaptivity. In Appendix B we present a graphical overview of the solution procedure for the modified Helmholtz equation.

Formulation
Consider the forced isotropic heat equation subject to initial-and Dirichlet boundary data U 0 and g, respectively. To fix notation let Ω be a time independent, compact (N Γ + 1)-ply connected region in R 2 with a boundary Γ consisting of (N Γ + 1) closed curves. These form the set Γ = {Γ n } N Γ n=0 , where Γ 0 is the outer boundary of the region Ω, see Fig. 1. The component curves are individually smooth and parametrisation each is assumed to be known. The outward directed normal at y ∈ Γ is denoted ν(y) = ν y and κ(y) denotes the curvature at y ∈ Γ.

Discretising in time and the modified Helmholtz equation
The heat equation (1) is first discretised in time, an approach known as elliptic marching or Rothe's method. To prevent severe time step restrictions an implicit-explicit (IMEX) scheme is used. It consists of using an implicit discretisation of the stiff terms and an explicit one for the nonstiff terms [30]. Regardless of the specifics of the IMEX scheme, to advance the solution U in discrete time a sequence of modified Helmholtz equations are solved. The modified Helmholtz equation is stated as with u unknown in Ω. The scalar parameter α 2 is inversely proportional to the time step δt; its explicit form along with f and g are given by the specific IMEX scheme. We use an adaptive IMEX Runge-Kutta method of fourth order in this paper, see Appendix A.1. However, what follows holds for any IMEX scheme.
Using the linearity of the differential operator α 2 − ∆, the solution u to (4)-(5) is decomposed into a homogeneous solution u H and a particular solution u P , such that u(x) = u H (x) + u P (x) for x ∈ Ω. First the particular solution is acquired by solving a free space problem assuming the existence of an extension f e ∈ C k (R 2 ), for some k ≥ 0, of the right hand side f from (4), such that for some finite L. The boundary condition, given by the Dirichlet data g in (5), is satisfied by u if u H is a solution to In short, first solve the free space problem (6) to obtain the boundary data for the homogeneous problem (10)- (11). The solution to the modified Helmholtz equation is the sum of the two solutions, u(x) = u H (x) + u P (x) for x ∈ Ω. See the flowchart in Appendix B for a graphical overview.

The inhomogeneous modified Helmholtz equation
The free-space modified Helmholtz equation (6)-(7) can be solved with Fourier transforms. Letû P =û P (ξ) and f e =f e (ξ) denote the Fourier transforms for u P and f e , respectively. Here ξ = [ξ 1 , ξ 2 ] ∈ R 2 with ξ = |ξ|. Then under the Fourier transform (6) is and we obtainû P (ξ) =f e (ξ) Note that the above expression is free of singularities, since α 2 0. The solution is given by the inverse Fourier transform For this solution to be well-defined the extension f e must be in L 1 (R 2 ). How to construct said extension and compute an approximation of u P is presented in 3.1.

The homogeneous modified Helmholtz equation
Consider the homogeneous modified Helmholtz equation (10)- (11). The free-space Green's function G(x, y) for the operator α 2 − ∆ is where K 0 denotes the zeroth-order modified Bessel function of the second kind. In other contexts the kernel G(x, y) is also referred to as the Yukawa or screened Coulomb potential. As in [1,4], we seek the solution u H (x) for x ∈ Ω in the form of a double layer potential: with the kernel where K 1 denotes the first-order modified Bessel function of the second kind. The limiting value as x goes to y along a boundary segment Γ n is well defined: where κ(y) is the curvature of Γ n at y ∈ Γ n , n = 1, . . . , N Γ . The density µ : Γ → R is not known a priori; it is found through the solution of a boundary integral equation. Such an equation of the second kind for µ can be formulated as For a derivation see e.g. [31]. Forg ≡ 0 only the trivial solution µ ≡ 0 along Γ satisfies (19). Thus by the Fredholm alternative the solution µ exists and is unique for any integrableg, for both simply and multiply connected domains [32]. This property is inherited by the corresponding discretised systems as well, introduced in Section 3.2. Each contour Γ n is split into N P,n intervals, referred to as panels, where Γ n,k is the kth panel on the nth contour and N P the total number of panels over Γ. A panel Γ n,k is represented by a known parametrisation γ n,k , such that By introducing a speed function s n,k (t) = |γ n,k (t)| and µ n,k (t) = µ(γ n,k (t)) the layer potential (16) can be written as and analogously for the boundary integral equation (19) µ

Discretisation
This section covers the numerical treatment of the modified Helmholtz equation. Note that two different methods are needed, one for the inhomogeneous problem and one for the homogeneous problem. We assume some suitable 5 IMEX scheme has been chosen for temporal discretisation of the heat equation (1)-(3), e.g. the Runge-Kutta methods presented in Appendix A.1.
Consider a box B = [−L, L] 2 in R 2 that containsΩ. The complement ofΩ relative to B is denoted by E. Denote the grid by X, which is a set of N 2 u elements x, referred to as nodes or points. They are uniformly distributed with spacing δx over B. Let subscripts indicate subsets of X, such as X Ω = {x ∈ X|x ∈ Ω} and X E = {x ∈ X|x ∈ E}.
X E = {x ∈ X | x ∈ E} such that it satisfies (8)- (9). Thereafter we consider the homogeneous problem (10)-(11), formulated as a boundary integral equation on Γ. The solution is computed at the locations X Ω = {x ∈ X | x ∈ Ω} in a post-processing step. The solution to the modified Helmholtz equation is computed at all grid points that fall inside Ω, i.e. the elements of X Ω . First we present how to find this solution for the free space problem (6)- (7). This involves extending the function f , based on the data at X Ω = {x ∈ X | x ∈ Ω} to

The inhomogeneous problem and function extension
An approximate solution to the free-space problem (6)- (7) is computed by discretising the integral in (14) with the trapezoidal rule. It is evaluated efficiently with FFTs on the regular grid X in B, thus in X Ω as well, and on the boundary Γ with a non-uniform inverse FFT. The latter is used to modify the given Dirichlet boundary data (5) for the homogeneous modified Helmholtz equation.
If the compactly supported f e in (6) is smooth, then the coefficients in the Fourier series expansion decay exponentially fast with the wave number, and this procedure would be specially accurate. With limited regularity, the Fourier coefficients instead decay algebraically, with one additional order for each continuous derivative. This approach requires an extension f e of f defined on X, preferably with high global regularity and compact support. It is constructed with PUX, which is briefly reviewed in this subsection. The basic concept is to blend local extensions by a partition of unity into a global extension with compact support, enforced by weight functions. The global regularity of the extension is directly related to the construction of said partition of unity. This is achieved by distributing overlapping partitions along the boundary Γ of Ω. In each partition the local values of f are used to extend it to the points in the partition that fall outside Ω. For a more extensive treatment see the original work [9].

Partition of unity
Let . The superscript k indicates the smallest subset C k 0 of C 0 that ψ k is a member of. Define a partition Ω i as the support of ψ k i , i.e. Ω i = supp(ψ k i ), which is a disc with radius R. We will return to the choice of ψ k in Section 3.1.3. Note that all partitions have the same radius. The number of partitions N ψ , the location of the partition centres {p i } N ψ i=1 and radius R are chosen such that the partitions cover Γ and that the partitions overlap with approximately a radius. The following notation will be useful. Each partition Ω i has a set of points on the uniform grid within R of p i , which we denote X i , rather than X Ω i . It can be split into two disjoint subsets: Let N i denote the number of elements in X i . Analogously, let N i,Ω and N i,E denote the number of elements in X i,Ω and X i,E , respectively. See Figure 3 for a graphical example. Given a function f : Ω → R, the function values at the locations X i,Ω are used to create a local extension f e i . We will return to the construction of the local extensions in Section 3.1.2, but for now assume their existence.
For every partition Ω i and its associated radial basis function ψ k i define the corresponding weight function w i as which belongs to the space C k 0 . By construction the set of weights which is referred to in the literature as Shepard's method [33]. See Figure 2 for a visualisation. This construction is used to combine the local extension However, (25) is not used, as we want an extension that it is continuous or of higher regularity as it is extended by zero outside its support. Refer to the set of partitions as extension partitions and now introduce also the zero . They are included in the partition of unity definition (24) and distributed such that they overlap the extension partitions, but do not intersectΩ. The associated local extension f e i is set to be identically equal to zero for i = 1, . . . , N 0 ψ . Hence, as the zero partitions are blended with the local extensions in the first layer of partitions , the global extension will be forced to zero over the overlapping region. Therefore zero partitions should be placed such that f e has a controlled decay to zero and that the size of the overlap with extension partitions are about the same, see Figure 4. Thus the global extension will in these parts have the same regularity as w p , as given by the regularity of the compactly supported radial basis function ψ k . The extension f e of f is given by As ψ k we use one of the compactly supported Wu-functions, which are tabulated after their regularity k, see Table  1 or [34]. There are other options, but the Wu-functions have compact support and are simple to implement. Note that they have lower regularity at the origin, e.g. the Wu-function listed as C 4 is only C 2 at that point. Moreover, the k + 1th derivative of ψ k is of bounded variation. The partitions centres are set to be nodes on the regular grid that are the closets to be boundary, yet still in X i,Ω . Thus evaluation of weight functions at the origin is omitted and higher regularity is maintained. With this, we have described how local extensions are combined into a global one. It remains to construct the local extensions

Local extensions
We now return to the construction of the local extensions f e i for each extension partition i = 1, . . . , N ψ . The local extension f e i is created as a weighted sum of radial basis functions, that interpolates the values of f at x ∈ X i,Ω and is evaluated at x ∈ X i,E . The radial basis functions are denoted φ j (x) = φ( z j − x ), to distinguish them from the radial basis functions ψ k . The elements of the set Z = {z j } N φ j=1 ⊂ supp(ψ k i ) are the centres for the RBFs, whose distribution for now is left unspecified. The standard form of an RBF interpolant at a point x is where λ j are unknown coefficients to be determined. We use where ε is a shape parameter setting the width of the Gaussian. The smallest interpolation error is obtained when ε is small, yet nonzero, but no general value can be given [35]. With some abuse of notation let X i and Z refer to vectors with the members of respective set as elements. Then, following the outline of [36] (1 − r) 6 + (6 + 36r + 82r 2 + 72r 3 + 30r 4 + 5r 5 ) Tab. 1: Wu-functions ψ k ∈ C k 0 , with compact support in r ∈ (0, 1) [34]. Here (·) + = max (0, ·). The listed regularity excludes evaluation at the origin.
Consider a scenario when f is known for all nodes in X i , then the associated interpolation problem to (27) can be written as with f X i = f (X i ). If N i ≥ N φ then Λ can be solved for in a least-squares sense. However, this is an unstable problem for several reasons. First, the conditioning of the problem is heavily dependent on the shape parameter ε.
For small ε the interpolation weights Λ oscillate between positive and negative numbers of large magnitude [35]. Furthermore, it is not uncommon for the condition number for the interpolation matrix to be of order 10 18 or more. These characteristics are common for interpolation with radial basis functions. Additionally, the data is represented on a uniform grid; collocating at these locations is the worst possible setting for interpolation, as with polynomials. These shortcomings can be circumvented by avoiding collocation and considering a least squares problem instead. Note that all problems mentioned above are purely numerical artifacts. The function space spanned by Gaussians is indeed a good approximation space. Decouple the centres Z of the radial basis functions from X i and assume they are distributed in a near optimal way with respect to minimising the interpolation error. We wish to omit explicit use of the interpolation coefficients Λ in (29). It can be achieved by formally solving for Λ by collocating at the centres Z: Here f Z are the values of f at the locations Z, which are unknown. Due to the choice (28) the matrix Φ(Z, Z) is symmetric and positive definite, thus the inverse Φ(Z, Z) −1 is well-defined. We can now reformulate (29) as Henceforth we use the shorthand notation For the purpose of function extension, sort the data points in Ω i such that where the components are of length N i,Ω and N i,E , respectively. Consequently, A can also be rearranged and split into two block matrices Since f is known at X i,Ω it can replace the corresponding entries in f X i (31) with f i,Ω = f (X i,Ω ). For each partition we obtain the system, with f i,E = f (X i,E ) unknown. For each partition i the values f i,Ω are mapped to the nodes Z to obtain f Z . Thereafter we obtain f i,E , which is the local extension. That is: This approach allows us to use a non-uniform distribution of RBF centres which significantly improves the stability, but still lets the data be represented on the uniform grid. We also avoid explicit use of the interpolation weights Λ. It remains to address the notorious ill-conditioning of Φ, associated with the shape parameter ε set small. This is achieved by applying the algorithm RBF-QR. It is intended for a formulation as (34), since it computes A, rather than Φ −1 , which acts as a mapping of data from non-uniformly to uniformly distributed locations. Said algorithm performs a change of basis for A, and in process the condition number is reduced, see [37]. By the use of RBF-QR the restrictions of choosing ε is lifted.

Properties of PUX
Four parameters need to be set for the PUX algorithm: the shape parameter ε for the width of Gaussians (28) used as interpolation basis, the partition radius R, the length L for the computational domain B = [−L, L] 2 and N u , where N 2 u is the number of uniformly distributed nodes over B. The remaining parameters can be set based on these values. Here we give the most important relations. For a complete discussion see [9].
Due to RBF-QR the shape parameter can be set small without risk of suffering from ill-conditioning. A good value is ε = 2, but the error in solving the modified Helmholtz equation is relatively insensitive.
Let P be number of uniform grid points per partition radius, denoted as This measure is used to choose ψ k from Table 1, and the number N φ of basis functions (28) per partition. To see how P relates to ψ k , consider the convergence of the error in solving the modified Helmholtz equation (4)-(5), assuming that only resolving u P limits the accuracy. If f e is smooth then the error has asymptotically spectral convergence. However, the extension inherits the regularity of the weight function w. Recall that by construction w ∈ C k 0 for a fixed ψ k (23). Consequently the error has an asymptotic convergence of 4 + k, if the kth derivative of f e is of bounded variation. A Wu-function of high regularity is harder to resolve than one of lower regularity. This implies that given a resolution P the error in resolving the Wu-function may hamper the convergence. As in [9] we use the heuristic relation k = min for choosing ψ k . In Section 4 we confirm that (36) is a satisfactory estimate for an optimal ψ k given P. Creating a local extension involves solving the least-squares problem A i,Ω f Z = f i,Ω for f Z for some i. It should be sufficiently overdetermined in order to be a well-posed problem. Given a P the number of unknowns N φ should be set accordingly to obtain a certain ratio of knowns and unknowns. Still, P can be of such magnitude that N φ is larger than required to obtain good results and the least-squares problem is more stable and cheaper to solve if the unknowns are few. Thus if the available data is abundant it can be downsampled to reduce P, and therefore N φ . Let c be the sampling parameter, defined as If c = 1 then all points are used, c = 2 means that every other point is removed, etc. Then, as in [9], we use N φ = min 0.8π(P/c) 2 /4, 3(P/c) (38) to set the number of radial basis functions per partition. Note that choosing ψ k is question about resolution; Wufunctions of higher regularity require larger P to be well resolved, while setting N φ is related to solving a least-squares problem. These are two separate problems and two different values for P may be used. So given a P we set ψ k according to (36) and then compute c with (37). Now N φ is set by (38) for P/c. Thus the local least-squares problems are solved on a potentially coarser grid, but the local extensions are on the original grid. The distribution of RBF-centres Z can be chosen freely, and we use the quasi uniform Vogel node distribution defined as in a unit disc. See Figure 3 for a visualisation. The distribution (39) is near optimal and RBF-QR performs well up to about 400 nodes. The locality of the weight functions guarantees that the least squares systems are of moderate size, which can be solved in parallel.
Constructing A (33) with RBF-QR is a computationally expensive operation, so employing it for every partition is undesirable. However, the matrix is the same for all partitions since p i is centred at a grid point from the uniform distribution. Thus the pairwise distances for the elements in X i are independent of i. Therefore a single matrix A can be precomputed with RBF-QR and reused for all extension partitions. The only difference between them in terms of A is the decomposition of X i into X i,Ω and X i,E , as it depends on how the boundary Γ intersects the partition. Note that the zero partitions may individually have a radius different from R in order to conform to the geometry of Ω and to overlap the extension partitions properly.

The homogeneous problem
For simplicity, assume the number of countours N Γ to be one and write Γ n,k = Γ k , y n,k = y k and s n,k = s k . We apply an N Q -point, panel-based Nyström discretisation scheme based on the composite Gauss-Legendre quadrature rule, with nodes t G m and weights W G m , with m = 1, . . . , N Q . Let y k,m = γ k (t G m ), s k,m = s k (t G m ) and µ k,m = µ k (t G m ). An approximation of the solution µ to (22) and correspondingly for (21) we have An important observation is that the kernel M (17) is not smooth and can contain singularities, depending on how x approaches y. Here the Gauss-Legendre quadrature rule is insufficient, as the resulting loss of accuracy can be critical enough to render the result useless. We elaborate on this topic in Section (3.2.1).
In matrix notation (40) can be written as (I + M)µ =g, where I is the identity matrix and M a compact operator. The density µ can be efficiently obtained with GMRES, in terms of numbers of iterations. The condition number for I + M is typically small or moderate and uniformly bounded. A fast multipole method (FMM) can be used for efficient computation of the involved potentials in (40) and (41) [38]. We use the point to point FMM for the two-dimensional Yukawa kernel presented [4]. It is based on the volume equivalent in [39]. For the corresponding three-dimensional version see [40].
Finally a note on the restriction of the boundaries being smooth. For non-smooth boundaries the integrand of (22) is not compact and the Fredholm alternative fails. While there are theoretical results on the solvability with Lipschitz continuous boundaries [41], they require the implementation of sophisticated quadrature techniques, such as [42], which we have not implemented. These methods also allow cusps, i.e. non-Lipschitz boundaries, and mixed boundary conditions.

Special purpose quadrature
When solving for µ in (40) or evaluating the layer potential (41) several orders of accuracy may be lost, since the kernel M (17) is not smooth. Moreover, M can be singular, depending on if x approaches some y ∈ Γ along Γ or from Ω. One of the most efficient methods to circumvent this loss of accuracy is explicit kernel-split quadrature with product integration by Helsing, see [19]. However, for the modified Helmholtz equation with large α, i.e. for high temporal resolution, it can fail completely. Below we sketch the problem, its relation to α and how to circumvent it.
We start by explaining product integration, which requires the involved integrals to be expressed in complex notation. To keep these paragraphs brief and simple, the reformulations are omitted. Consider a single panel Γ k ∈ C with endpoints at −1 and 1, but the panel does not have to follow the real axis. Let ϕ : Γ k → R be a smooth function and s : Γ k × C → R a non-smooth kernel that may be singular or nearly singular. The goal is to compute accurately for some fixed τ 0 ∈ C arbitrarily close to or on Γ k . To do this, approximate ϕ with a polynomial of degree with unknown coefficients {c n }. Inserting this into (42) gives The integrals on the right hand side can can be computed analytically through recursive formulas. The unknown coefficients {c n } are obtained by solving a Vandermonde system. If ϕ can be accurately represented as a N Q − 1 degree polynomial over Γ k , then product integration allows evaluation of integrals such as (44) without loss of accuracy as τ 0 and τ approach each other. Kernel-split means that a kernel is decomposed into smooth and singular terms. Leaving complex notation, by [43, §10] the first-order modified Bessel function of the second kind K 1 , appearing in (17), can be decomposed as This form is attractive since the singular terms are separated and can be studied individually. Here I 1 is the modified Bessel function of the first kind of order one and K S 1 is a power series in x. For the kernel M, see (17), the situation is slightly more involved, as the singularity structure depends on how x approaches y ∈ Γ. To distinguish between the two cases, for any y ∈ Γ denote M(x, y) as M Γ (x, y) for x ∈ Γ and M Ω (x, y) for x ∈ Ω. We first study M Γ ; the decomposition (45) motivates the formulation with M Γ,L identified as The term M Γ,0 is smooth and by (18) we have since the term log( x − y )M Γ,L (x, y) goes to zero in the limit x → y. But in this limit the derivative of log( y − x )M Γ,L (x, y) has a log-type singularity. Thus standard quadrature rules that relies on smoothness fail to be accurate.
To maintain accuracy product integration is needed, even though the limit is well-defined. In terms of (44) φ and s correspond to µM Γ,L and log. This approach is used to compute the involved integrals in (22).
In the case x ∈ Ω, corresponding to computing (21), the kernel M(x, y) is singular in the limit x → y and product integration is required. We have where M Ω,0 is a smooth function, M Ω,L = −α 2 /2M Γ,L and M Ω,C = −α 2 /2. Again, we identify ϕ from (44) as µ multiplied with M Ω,L or M Ω,C and the singular function s corresponds to either log( y − x ) or (y − x) · ν y / y − x 2 . In complex notation, the latter is reduced to a Cauchy-type singularity. Both M Γ,L and M Ω,L contain the factor I 1 (α x − y ), which grows like e α x−y / α x − y . The scaling with α can make I 1 grow too fast over a single panel to be accurately approximated by e.g. a 15th degree polynomial or even a 31th degree polynomial. The product integration relies on φ being well approximated by such a polynomial (43), otherwise the result may be very inaccurate. An adaptive time stepper will adjust the time step to satisfy the given tolerance, potentially decreasing it until the algorithm stalls.
This problem is not unique to the modified Helmholtz equation, but appears for biharmonic and Stokes equations as well. One solution is an algorithm presented in a separate paper, see [18]. By local refinement of panels through adaptive recursive bisection a kernel-split quadrature with product integration can be used successfully for a wide range of α. It is ensured that the new panels are of adequate size to accurately approximate φ with polynomial interpolation. The method is effective in terms of computations, as the increased cost scales as log(α). Moreover, K 1 (α x − y ) ∼ π/(α x − y ) e −α x−y for large arguments, i.e. K 1 is very localised for large α and only a small portion of the boundary Γ needs to be upsampled.

Numerical results
In this section we present numerical results, starting with a study of the modified Helmholtz equation to confirm that the parameters for PUX can be set as in [9] for the Poisson equation. It forms the basis for the second numerical experiment, where the modified Helmholtz equation is solved on a more complex domain. The heat equation is solved on the same domain, for a range of set tolerances with an adaptive time stepper for different grid resolutions. Finally, the Allen-Cahn equation , a reaction-diffusion problem, is solved with randomised initial data.
To compute the errors we consider an evaluation grid. It consists of N 2 eval uniformly distributed nodes over the computational domain B. We evaluate the numerical solution and an analytical or computed reference solution on the nodes that fall inside Ω. The cardinality of this set of nodes as N eval,Ω . Two different errors are computed: the relative and for a vector u of length N eval,Ω . When referring to the errors we mean both of them.
The following parameters are user specified in the numerical experiments: the length L for the computation domain B = [−L, L] 2 , the resolution N u , partition radius R and the number N P,n of Gauss-Legendre panels for each component curve Γ n . We use set the shape parameter ε = 2 for all numerical experiments and set the number of Gauss-Legendre nodes N Q = 16.

Example 1: Study of weight functions
We now solve the modified Helmholtz equation (4)-(5) for u(x, y) = sin(2πx) sin(2πy) exp(−(x 2 + y 2 )), to confirm that the parameters N φ and c and the function ψ k can be set by (35), (38) and (37), as in [9] for the Poisson equation. To reduce the complexity of the problem, assume the corresponding right hand side to be known in all of R 2 , not just Ω. To isolate the influence of the choice of weight function ψ, see Table 1, the actual values of f are used as values for the local extensions f i,E , instead of the extrapolated ones A i,E f X i . Compact support is still enforced via PUX, but blending with the zero partitions reduces the regularity of f e to k. The computational domain is the unit circle centred at (17/701, 5/439), contained in the box B = [−L, L] 2 , with L = 1.5. The resolution N u attains values between 40 and 500 and for the evaluation grid use N e = 1000. The partition radius and the number of panels are set such that only the resolution of the uniform grid X limits the accuracy. In this case the partition radius is R = 0.4 and the number of panels N P = 32. This means that the rate of convergence is only dependent on the regularity of the extension and we can study the influence of choice of Wu-function. Furthermore, we set α 2 = 10.
In Figure 5 the errors for solving the modified Helmholtz equation are plotted as functions of the number of grid points for different Wu-functions. The behaviour of the errors is as for the Poisson equation in [9]: ψ k with few continuous derivatives requires less points to be represented then ψ k with a larger k. Consequently, high regularity can increase the error, since ψ k is not sufficiently resolved. Compare the errors for using ψ 1 and ψ 5 in Figure 5 for N u ∼ 40. As the grid is refined the decay is spectral until the errors is limited by an algebraic tail. The algebraic tail has a slope of 4 + k, as expected. The ∞ -error is about one to two digits less accurate than the 2 -error, which is consistent for all numerical experiments in this paper. The reason is that there is almost always some target points close to the boundary for which the special quadrature does not give optimal results, e.g. at the intersection of two panels.
We now solve the modified Helmholtz equation in the same numerical setting, but let ψ k be set automatically by (36). The result is presented in Figure 6 and the lines follows the corresponding lowest errors in Figure 5. Thus (36) indeed chooses ψ k correctly for a given N u and we can set the PUX parameters for the modified Helmholtz equation as for the Poisson equation. This holds for α 2 from 10 to 10 5 as well, as is shown in the following numerical experiment. Moreover, the error decreases as that of a tenth order method. For the subsequent numerical experiments ψ k , N φ and c are set by (36), (38) and (37)

Example 2: the modified Helmholtz equation on a multiply connected domain
We now study the modified Helmholtz equation with a more complex setup for α 2 = 10 n , n = 1, 2, 3, 4, 5. We take the solution to be u(x, y) = cos 20 on the multiply connected domain shown in Figure 7 and evaluate the right hand side in (4) accordingly. The corresponding extension by PUX is shown in 8, where N u = 1000 and k = 5. The outer boundary is discretised into 80 panels, and the boundaries of the cavities are discretised with 20 panels each. Again all parameters are set such that only N u sets the bound for the error. We set R = 0.23 and L = 1.2. The parameters ψ k , N φ and c are set by (36), (38) and (37). The parameter α 2 ranges from 10 to 10 5 . The evaluation grid has a resolution of N eval = 1000. The results in Figure 9 suggest that (36) is a good estimate for setting ψ k for more complex problems as well. We obtain 10th order convergence with grid refinement. Note that slightly better results can be achieved; the same parameters are used for the entire range of α and are therefore potentially not optimal. As in the previous example the relative ∞ -error is about two orders of magnitude larger than the relative 2 -error. In Figure 7 the largest error is by the rightmost point in Ω, at the intersection of two panels. The special purpose quadrature is know to struggle with maintaining full accuracy in such situations.
The modified Helmholtz equation becomes significantly harder to solve for increasing α 2 . This is due to the rapid decay of the kernel (17), which requires a very fine resolution of the boundary to be resolved. We also suffer from cancellation errors due to the scaling of terms with α or α −1 . Still, this is not alarming, as an relative ∞ -error of about 10 −10 can still be obtained for α 2 = 10 5 . In terms of the heat equation this corresponds to a time step of about 10 −5 .

Example 3: Adaptive time stepper
We now test the solver for the heat equation (1)-(3) by setting a tolerance for the time stepping error and investigate if it can be maintained for different resolutions N u . For this purpose we use the IMEXRK34 scheme with an adaptive time stepper, see Appendix A.1.3 and Appendix A.1.2. The smaller the time step, the harder the modified Helmholtz equation is to solve, as concluded above. Thus a high order time stepping scheme, such as IMEXRK34 of fourth order, is a suitable choice since larger time steps can be used. However, other time marching methods can be used as well.
The domain and all parameters are set as for the previous experiment. The heat equation (1)- (3) is solved with right hand side F, initial condition and Dirichlet boundary data prescribed by the analytical solution U(t, x, y) = exp(−t) sin((x cos(π/4) + y sin(π/4))) + cos 20 where the time ranges from 0 to 1. For the evaluation grid we set N eval = N u and measure the error at terminal time t = 1. In Figure 10 the red lines correspond to set tolerances. It is clear that the adaptive time stepper works as intended, even for tolerances down to 10 −10 . The relative ∞ -error is more sensitive to the resolution and exceeds the set tolerance earlier in terms of spatial resolution, roughly with one order in magnitude.

Example 4: The Allen-Cahn equation, a reaction diffusion problem
The Allen-Cahn equation is stated as with C = 10 −3 . The right hand side of (55) is nonlinear and has three stationary points: U = −1, 0, 1. For randomised initial data the solution creates over time patterns with zones attaining these values. The initial data is not entirely randomised, since we need smoothness to discuss convergence and accuracy. Instead, we create smooth data by uniformly distributing 50 Gaussians (28) with ε = 10 over the computational domain with L = 1.2. Each Gaussian is assigned a coefficient drawn randomly from a uniform distribution over −0.5 to 0.5. The partition size R is set to 0.1; the domain, the extended right hand and the distribution of partitions are shown in Figure 11. Each boundary component is discretised with 80 panels. We create a reference solution by solving the Allen-Cahn equation with tolerance 10 −6 and N u = 800, from time 0 to 6. The errors are measured on grids with N eval = 200, 400 at terminal time t = 6. Snapshots of this solution are shown in Figures 12a to 12f. Indeed the solution forms a pattern of patches with the values −1, 0 and 1. The results are shown in Table 2. For N u = 400 the relative 2 -error stays under the set tolerance. However, unlike example 3 the relative 2 -error is always a factor ten larger than the set tolerance. For N u = 200 only the tolerance 10 −3 can be obtained. Clearly this resolution is insufficient to resolve the spatial problem more accurately than that. The error at the terminal time t = 6 for N u = 400 with tolerance 10 −5 is shown in Figure 13. In this figure the evolution of the time step is also shown; as the solution advances in time the time step becomes larger. Initially it grows faster, compared to later, as the initial time step was intentionally set small. The initial data U 0 (56). Right: The right hand side of (55) at t 0 , extended with PUX. Black corresponds to zero partitions and red to interpolation partitions. Note that to increase visibility of the field a different scaling is used than for 12a-12f

Conclusions
We present a framework built around a panel-based Nyström boundary integral method for solving the forced isotropic heat equation in two dimensions, on multiply connected complex domains. We have addressed several of the issues listed in [1], thereby increasing the class of solvable problems as well as the accuracy in the solutions.
We show how any IMEX method can be applied as time stepping scheme, and employ an adaptive fourth order Runge-Kutta scheme in our examples, to accurately solve the heat equation as well as the Allen-Cahn equation, a reaction-diffusion problem with a nonlinear forcing term. Regardless of the specific details of the chosen method, a time step in solving the heat equation is reduced to solving one or a sequence, for a multi-stage method, of modified Helmholtz equations.
As in [4] we formulate the modified Helmholtz equation as a boundary integral problem. Utilising the linearity of the differential operator, the solution is split into a particular-and homogeneous problem. Solving the former to high accuracy relies on extending the given right hand side from the domain it is given on to the entire plane. It is achieved with a partition of unity extension (PUX), that only requires known data at uniform point locations inside the domain. The extension that is computed on a uniform grid in a rectangular domain has compact support and a specified global regularity, making spectral methods very efficient and simple to use. We confirm that the various parameters for PUX, in the context of the modified Helmholtz equation, indeed can be set as for the Poisson equation in [9]. This yields an automated selection for the global regularity to balance different errors, leading to a method which converges with an order 10 in the grid size.
A panel-based Nyström boundary integral method is used to solve the homogeneous problem with modified Dirichlet data, such that the total solution is the sum of the particular-and homogeneous solution. The boundary values of the particular solution are computed using a non-uniform FFT. For evaluation of singular and nearly singlar integrals, we have introduced a methodology based on product integration and an explicit kernel split that has given highly accurate results for the Helmholtz [15] and Stokes equations [17]. For large α (small time steps), the method in its original form, would however fail completely if an unfeasibly high upsampling of the boundary was not applied. We however realized that this upsampling is only needed very locally, and developed an adaptive approach [18] to achive a computationally efficient method with high accuarcy.
In total, these developments yields a method for very accurately solving the heat equation on comlex domains. The highest attainable accuracy in the solution of the modified Helmholtz equation does show a weak depenence on α, but even for the largest values, solutions can typically be attained with at least ten correct digits, meaning that strict time stepping tolerances for the heat equation can be satisfied.
In terms of future developments, it would be useful for some problems to replace the uniform grids and FFT-based method for the particular solution with a volume potential evaluation based on an adaptive FMM. This would however need an integration of the PUX method into the adaptive procedure. Another development is to consider the solution of the heat equation and the closely related advection-diffusion equation on time-dependent domains. The motivation for this is the need to solve such an equation for the concentration of surfactants in the oil-phase of a micro-system with water drops in oil. These surfactants, or surface active agents, have an exchange with surfactants on the drop surfaces, that alters the surface tension of the drop. Numerical methods for simulating surfactant advection and diffusion on the boundary of drops have been understood and implemented successfully, see [44,45]. An important extension would be to allow also for surfactants in the oil-phase. One strength of these methods is the accurate treatment of interface conditions, something that is absolutely essential at these small scales where the interface dynamics is of key importance.

Acknowledgements
We thankfully acknowledge the support of the Swedish Research Council under Grant No. 2015-04998 and funding from the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine. We are humbly grateful for the support from the Natural Science and Engineering Research Council of Canada.
Appendix A. Adaptive time-stepping with IMEX Runge-Kutta methods Appendix A. 1

. Adaptive discretisation in time
This appendix shows how applying implicit-explicit Runge-Kutta (IMEXRK) schemes from [5] to the heat equation reduces it to a sequence of modified Helmhotlz equations to solve at each time step. Formulate the heat equation (1)-(3) as where the superscripts denote implicit and explicit, referring to the term being classified as stiff or nonstiff, respectively. Let t N denote an instance in time that is the sum of previous discrete time steps {δt i } N i=1 that may be of different size: for some initial time t 0 . Let U N be the approximation of U(t N ), then the approximated solution at time t N+1 is where N S is the number of stages for k σ , σ ∈ {I, E}, computed as The second argument of F σ in (A.5) is defined as a σ i, j k σ j + δt N+1 a I i,i k I i , i > 1, (A.6) andŪ 1 = U N . The coefficients {a σ i, j } N S i, j=1 , {b σ j } N S j=1 and {c σ i } N S i=1 are tabulated in the two associated Butcher tableaus for σ = I and σ = E, see Table A.3 for a general IMEXRK scheme. The principal difference between the coefficients for implicit and explicit methods is that a E i, j = 0 for i ≤ j while a I i, j 0 for i = j, excluding i = 1. The quantityŪ i is unknown for every i = 2, . . . , N S , since the corresponding implicit stage k I i is unknown. The implicit stage at i is k I i = F I = ∆Ū i by definition (A.2). To avoid approximating the differential operator replace k I i in (A.6) with ∆Ū i and reformulate as The idea is to solve forŪ i and since the right hand side is known, ∆Ū i can be extracted from the expression above. The equation (A.7) has the form of the modified Helmholtz equation (4)-(5): f (x) corresponds to the right hand side, u(x) =Ū i (x) and α 2 = (δt N+1 a I i,i ) −1 . We stress that α 2 ∼ (δt N+1 ) −1 ; the larger α 2 is the harder (4)-(5) is to solve accurately in terms of numerics, see Section 3.2.1. The associated boundary condition g is (3) evaluated at t N +δt N+1 c I i . To obtain the next stage k  · · · · · · a σ N S N S b σ 1 · · · · · · b σ N S Tab. A.3: Coefficients for an IMEXRK scheme, where σ ∈ {I, E}, denoting implicit or explicit, applied to the stiff and nonstiff term, respectively. In general a E i j = 0 for i ≤ j and a I i j 0 for i = j, excluding i = 1.
With k I i known the stage k E i , that is F E , can be computed explicitly. Note that for (1)-(3) F E = F(t, x), so the explicit stage k E i is independent of the implicit stages, thus it is computed directly. Note that this is not the case if e.g. an advection term ∇U is added, as it would be included in F E . In order to keep the formulation general we think of F E as function of U.
To summarise: the approximate solution U N+1 at time t N+1 is given by (A.4). The stages k I i , for i = 1, . . . , N S are obtained by solving (4)-(5), corresponding to (A.7), and explicit computation of (A.8). OnceŪ i is known k E i = F E t N + δt N+1 c E i ,Ū i is computed explicitly. See the flowchart in Appendix B for a graphical overview.
Given U N , the next solution U N+1 is obtained by (A.4). The implicit stages (A.5) must be solved for. Explicit stages are computed directly. The same stages are used to compute low order approximationÛ N+1 , used for adaptive time stepper.
Input: time step δt N+1 , Butcher tableau A.3, solution U N at t N , first implicit stage k I 1 and Dirichlet boundary data.
Solve for U I i from (A.7), by solving (4) Are all stages computed?
Is r < TOL?
Output: Solution U N+1 , time t N+1 = t N + δt N+1 , time step δt N+1 and k I 6 , which is k I i for the next iteration in time.
Approximate solution U N+1 to the diffusion equation at time t N+1 . Solve modified Helmholtz equation u(x) = g(x), x ∈ Γ.
Decompose u = u P + u H .
Input: α 2 , f and g Construct extension f e (x) of f with PUX, see Section 3.1.
Solve α 2 u P − ∆u = f e in Fourier space with FFT, Section 3.1. Compute u P in Ω and u P | Γ on Γ.
Solve α 2 u H − ∆u H = 0 in Ω with u H =g on Γ as in Section 3.2.
Compute u H in Ω The solution to the modified Helmholtz equation is u = u P + u H .
Output: u Donẽ g = g − u P | Γ u H u P