{'fraction_non_alphanumeric': 0.04067866964146734, 'fraction_numerical': 0.014930213251699269, 'mean_word_length': 4.147426053281698, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An integral equation method for the advection-diffusion equation on time-dependent domains in the plane

Boundary integral methods are attractive for solving homogeneous linear constant coefficient elliptic partial differential equations on complex geometries, since they can offer accurate solutions with a computational cost that is linear or close to linear in the number of discretization points on the boundary of the domain. However, these numerical methods are not straightforward to apply to time-dependent equations, which often arise in science and engineering. We address this problem with an integral equation-based solver for the advection-diffusion equation on moving and deforming geometries in two space dimensions. In this method, an adaptive high-order accurate time-stepping scheme based on semi-implicit spectral deferred correction is applied. One time-step then involves solving a sequence of non-homogeneous modified Helmholtz equations, a method known as elliptic marching. Our solution methodology utilizes several recently developed methods, including special purpose quadrature, a function extension technique and a spectral Ewald method for the modified Helmholtz kernel. Special care is also taken to handle the time-dependent geometries. The numerical method is tested through several numerical examples to demonstrate robustness, flexibility and accuracy


Introduction
Many physical phenomena of interest in science and engineering involve moving domains. Yet there is no standardized numerical solution method for such problems, thus solving partial differential equations (PDEs) on time-dependent geometries is an active area of research. The main problem for classical numerical methods is that the underlying grid, where the data is represented, must be remeshed to follow the movement of the domain. Clearly, this becomes expensive as the temporal and the spatial resolutions are increased. Instead, it is desirable to keep the underlying grid fixed, which comes with its own set of challenges, depending on which method is used. We will begin by introducing the problem and our proposed method, and with this as background we will then discuss other methods.
We propose a boundary integral-based method for solving the isotropic advection-diffusion equation, also called the convection-diffusion equation, on time-dependent geometries in two spatial dimensions. The advection-diffusion equation is formulated as where U is an unknown quantity to be solved for, v is a velocity field that is either given explicitly or is the solution to an additional PDE, the constant D > 0 is the diffusion coefficient, the source term J(t, x) is a smooth explicitly given function, Γ(t) is the boundary of the domain Ω(t) ⊂ R 2 at time t with at least two continuous derivatives,n is the outward directed unit normal at x ∈ Γ(t) as shown in figure 1, g is given Neumann boundary data, and the domain Ω(t) = Ω(t)∪Γ(t) is bounded and continuous in t. We assume that the velocity field v transports U and Γ with the same velocity. Multiply connected domains are not addressed in this paper, but including them requires no additional theory or computational methods. The extension to exterior problems is straight forward, assuming suitable far-field conditions are supplied. We do however consider the advection-diffusion equation above with periodic boundary conditions for parts of the boundary. The advection-diffusion equation is often part of models of complex problems in computational physics, for example the dynamics of multiphase flows on the microscale with soluble surfactants (surface-active agents). Surfactants decrease the surface tension at interfaces between fluids. At the microscale these forces dominate, as the surface-to-volume ratio is large, thus surfactants have a significant impact on the dynamics of the multiphase flow. They can be used to control droplet dynamics, which is central for many microfluidic applications [51]. Surfactants are usually present on the interface and in one of the fluid phases. This can be modeled with the advection-diffusion equation, where boundary interface conditions describe the exchange between the drop interface and the bulk. Considerable effort is given to accurately simulate microfluidic flows, as computer simulations is a much used tool to study such problems.
Here boundary integral methods enter the picture. The velocity field in drop and particle suspension problems at the microscale is often modeled with the Stokes equations, which can be recast as a second-kind Fredholm integral equation [44]. Consequently, the dimensionality is reduced as the unknowns are confined to the boundary of the domain, including particle and drop surfaces, the latter being fluid-fluid interfaces. Thus only the boundary needs to be discretized [46,50]. Tracking the evolution of the interfaces and the interface conditions are naturally handled with a boundary integral method, effectively avoiding remeshing the domain as it evolves in time. Another advantage of a boundary integral formulation is that high accuracy can be maintained as droplets, vesicles or particles get arbitrarily close to each other.
Time-dependent PDEs can be approximated by a sequence of inhomogeneous elliptic PDEs under appropriate temporal discretization, such as implicit-explicit (IMEX) spectral deferred correction (SDC) or IMEX Runge-Kutta methods (R-K) [11,4]. The resulting PDEs can be solved with a boundary integral method, which has been addressed for stationary geometry [18,3]. In the case of the advection-diffusion equation the diffusion term is treated implicitly, reducing the equation to a sequence of modified Helmholtz equations, ∂u(x) ∂n = g(x) , x ∈ Γ, (6) where u corresponds to an approximation of U at the next instance in time, α 2 = O 1 δt with time-step δt, and the right-hand side f is updated at every time-step. We write the solution to (5), (6) as a sum of a single layer potential and a volume potential, the latter is due to the PDE being inhomogeneous for f = 0.
Computing the required volume potential poses several numerical challenges if Ω is not geometrically simple, e.g. rectangular. We take the approach of creating a continuation f e of f , with a specified regularity and compact support on a box B embeddingΩ, such that f e | Ω = f . The extension f e is constructed with partition of unity extension (PUX) [17]. Now, let u P be the convolution of f e with the modified Helmholtz Green's function over B. The convolution can be computed efficiently with geometrically unaware fast summation methods, such as a fast Fourier transform (FFT). The solution to (5), (6) is the sum u P + u H , if u H satisfies the homogeneous counterpart of (5) with Neumann boundary data g(x) − ∂u P (x) /∂n. The solution u H is efficiently solved for with a boundary integral method, obtained via a Fredholm integral equation of the second kind, aided with special quadrature methods to maintain high accuracy for target points arbitrarily close to the boundary [19]. The decomposition of u into u P and u H allows us to enjoy the benefits of boundary integral methods for inhomogeneous elliptic PDEs.
When discretizing the boundary integral equation, and using an iterative method to solve the linear system, matrix-vector multiplications with a dense matrix must be rapidly computed. For the modified Helmholtz equation these matrix-vector multiplications are discrete sums with modified Bessel functions of the second kind of zeroth-or first-order. In a free-space setting, the sums can be computed directly at cost O(N 2 ) for N sources and targets. As the problem size grows, the evaluation becomes prohibitively expensive. There are several ways to speed up the computations, such as with a fast multipole method (FMM) [22,32] or a spectral Ewald method. The latter is an FFT based method, which is especially suitable for periodic problems. It has been implemented for the Laplace equation in 3D [1,40], the Stokes equation in 3D for both freespace and periodic setting [38,2], and in 2D for a periodic setting [45]. In this paper we present an Ewald method for the modified Helmholtz equation. An Ewald method is based on an Ewald decomposition, meaning that the sum is split into a "real space" sum which converges rapidly in R 2 , and a "k space" sum which converges rapidly in the frequency domain. By accelerating the evaluation of the latter using the FFT, the total computational cost of the method is O(N log N ).
When Ω evolves discretely in time the right hand side f in (5) may not exist everywhere in Ω. The function extension technique PUX is also used to address this problem of missing data. The basic concept of PUX is to blend local extensions by a partition of unity into a global extension. The regularity of the extension is directly related to the construction of said partition of unity. The local extensions are computed through interpolation with radial basis functions (RBF). The associated interpolation matrix is precomputed such that the ill-conditioning of the problem is significantly reduced.
The main contributions in this paper are the following.
• An accurate boundary integral method for the solution, and the gradient thereof, of the modified Helmholtz equation with Neumann boundary data over a wide range of values of the parameter α in (5).
• A spectral Ewald method for the kernel associated with the modified Helmholtz equation, both in a periodic and free-space setting.
• A high order solver of the advection-diffusion equation on moving domains in a boundary integral framework, that is adaptive in time.
The result is a solver for the advection-diffusion equation (1)-(4) of order K in time and tenth order in space, with K(K − 1) modified Helmholtz equations (5), (6) to solve per time-step.
An alternative boundary integral approach to the one presented in this paper is the one of Wang et al. [56]. Here, the solution is obtained through a direct approximation of the heat kernel, thus avoiding the discretization of the differential operator with respect to time. The corresponding layer potential contains no volume potential, but temporal integration from the initial time t 0 to some t. This integral is split into two: one "history part" from t 0 to t − δ, and one "local" part, integrating from t − δ to t. The former is computed with Fourier approximation and nonuniform FFT, which tends to be inefficient for small δ due to having to overresolve the frequency domain. The latter is computed with an FMM. The result is an efficient algorithm, that potentially can be extended to inhomogenoeus problems as well. A similar approach for the inhomogeneous counterpart, but for stationary domains, is found in [37], where the heat layer potential contains a volume potential. Here the layer potential is also composed of a local part and a history part In [15] Frachon and Zahedi present a Cut finite element method (CutFEM) for solving the Navier-Stokes equations on time-dependent geometry. In a CutFEM the geometry can be arbitrarily located on a fixed underlying mesh. The involved space-time integrals are in a variational formulation, and are computed by first using a quadrature rule in time. Then for each quadrature point the integrals in space are computed with CutFEM. In time, the trial and test functions are discontinuous between one time interval [t n , t n+1 ] to another. Continuity in time is enforced, weakly, by adding a correction term. Stabilization terms are also added which implicitly extend data to points in the new domain where data does not exist. Despite the addition of preconditioning and stabilizers, the method suffers from large condition numbers as the order of the quadrature scheme is increased. Thus while CutFEM methods are versatile in handling timedependent geometry, high-order in time is difficult to achieve. Another FEM-based approach is presented in [36] where an unfitted mesh is used in the spatial domain. A time-discretization based on finite difference approximation is used, rather than space-time variational formulation. Here extension methods are applied to address the problem of missing data when the domain is updated in time.
Yet another approach is a meshless method of high spatial order by V. Shankar et al. in [48]. Here the solution is represented by an approximation with radial basis functions, centered at a scattered node set. These nodes are used when applying finite differences. It then remains to solve for the interpolation weights to obtain the solution. The corresponding differentiation matrix does not need to be recomputed at each iteration in time, instead it is updated, which with their techniques is only required for nodes neighboring the boundary. Impressive results are presented for spatial convergence, but any high order methods for temporal accuracy are not presented.
The immersed boundary smooth extension method is used in [28] to solve the bulk advectiondiffusion equation on time-dependent geometry for simulating the Stefan problem. By applying an IMEX scheme they also obtain a sequence of modified Helmholtz equations to solve per timestep. Unlike the method presented in this paper, the function data f is not extended outside Ω, but instead an extension of the solution that satisfies the given PDE is solved for. This is achieved by imposing the unknown solution as boundary data for a homogeneous poly-harmonic PDE of degree k, for which the solution is the extension of regularity k. Third-order temporal and pointwise spatial convergence is demonstrated for a classical Stefan problem. This paper is organized as follows. First the mathematical formulation of the modified Helmholtz equation is presented in section 2, followed by section 3 covering time-stepping schemes, and numerical methods for the modified Helmholtz equation. Thereafter, PUX and its application to time-dependent geometry are introduced in section 4, followed by section 5 where numerical results are presented. Finally, conclusions and an outlook are given section 6.

Formulation
To solve the advection-diffusion equation (1)-(4) we apply elliptic marching, also know as Rothe's method or method of lines [27,9], meaning that we first discretize first in time, using an implicit treatment of the stiff diffusion term; the advective term and the source term are treated explicitly. Regardless of the specifics of the such discretization method, the result is a sequence of modified Helmholtz equations (5), (6) to be solved for each time-step. For example, a first-order IMEX Euler method applied to the advection-diffusion equation results in where . It can be formulated as the modified Helmholtz equation (5), . Note that f exists on Ω n , which may be different from Ω n+1 , and thus f is not defined in all of Ω n+1 . For the remainder of this section we will assume f to be defined both on Ω n and Ω n+1 . We will return to the question of how such an f is obtained in section 3.1.3. The free-space Green's function for the differential operator where K 0 is the zeroth-order modified Bessel function of the second kind, and · is the Euclidean norm unless specified otherwise. Note that neither the modified Helmholtz equation nor its Green's function are consistently named in the literature, but are also referred to as the screened Poisson equation, the linearized Poisson-Boltzmann equation, the Debey-Hückel equation and the Yukawa equation, and the Yukawa potential or screened Coulomb potential [29,8,47]. By the linearity of the differential operator α 2 − ∆, the problem (7) can be decomposed into a particular problem, with free-space boundary conditions u P (x) → 0 as x → ∞, and a homogeneous problem Thus one PDE includes the right-hand side f , also referred to as the forcing term, while the other enforces the Neumann boundary conditions with a correction ∂u P (x)/∂n. By construction U n+1 = u P + u H satisfies (7). The structure outlined above can be obtained using any IMEX rule. The methods differ in the definitions of α 2 and f , the required number of modified Helmholtz equations to solve per time-step, and how these results are combined to advance an approximate solution to the advection-diffusion equation in time. For example, with an IMEX R-K method the right-hand side f is a linear combination of intermediate so called stages. While this is the approach in [18] by Fryklund et al., where the heat equation is solved in a similar framework, here an IMEX SDC is used instead. The choice of SDC over R-K is discussed in section 3.1. The concept of SDC is to improve an initial approximation through solving a series of update equations with low-order methods, such as the first-order IMEX Euler method above. In the context of this paper, these update equations also take the form of the modified Helmholtz equation.

Particular solution
Here we present the mathematical formulation of the solution to the particular problem (9). Leaving numerics aside, a mathematically straightforward approach is to express the solution as a convolution of the forcing term f and the free-space Green's function (8) over Ω. Even though this formulation is simple, to numerically evaluate the integral accurately is complicated for complex geometries. Since we consider time-dependent geometries, high-order tailored quadrature rules for a specific geometry is not a feasible option. Therefore, we do as is proposed in [18,3], and explicitly write the solution as a Fourier series, after first reformulating the problem.
Instead of the initial formulation of the particular problem (9), consider where f e ∈ C q (B) for some q > 0,Ω ⊂ B and There is flexibility in choosing B. Thus, assuming the existence of an extension f e with compact support enclosingΩ, we may set B to be a box [−L/2, L/2] 2 containing supp (f e ), consider (11) as a periodic problem with periodic boundary conditions on B, and express the solution as a Fourier series. Standard results yield an explicit expression for the solution where {f e (k)} k , k = n2π/L with n ∈ Z 2 , are the Fourier coefficients in the Fourier series expansion of f e . Analogously, the gradient of the solution is Note that the zeroth-mode is well-defined for α = 0.

Homogeneous solution
Here we consider the homogeneous problem (10) in a boundary integral formulation which takes the geometry and the boundary data into account; the forcing term f is incorporated into the particular problem. The homogeneous solution in the form of a single layer potential is denoted by and the gradient of the solution can be written as The kernels are given by the Green's function (8), with the factor 1/2π extracted to outside the layer potential integrals, and where the gradient is with respect to x, and where K 0 and K 1 are the zeroth-and first-order modified Bessel functions of the second order, respectively [32]. Here, the layer density µ is assumed to be continuous and µ : Γ → R. It is the solution to the second-kind Fredholm integral equation with the double layer kernel The boundary integral equation (19) is obtained by enforcing the Neumann boundary condition in the limit as x ∈ Ω goes towards the boundary and using the jump conditions given in [32]. We note for future reference that the limiting value of D on the boundary is where κ is the curvature of Γ at y, assuming Γ ∈ C 2 . The kernel D is thus bounded and continuous along Γ, and therefore the integral operator is compact [30]. Since there is no nontrivial homogeneous solution, the Fredholm alternative states that (19) has a unique solution for any continuous data g(x) − ∂u P (x) /∂n [31].

Numerical methods
In this section we present the numerical methods used to solve the advection-diffusion equation (1)-(4). First we present the selected time-stepping scheme and how its application to the advection-diffusion equation results in a series of modified Helmholtz equations. Next, we demonstrate how to solve the modified Helmholtz equation on stationary domains. In section 4, moving domains are introduced, in combination with the PUX method for function extension.

Time-stepping
Here we describe an IMEX SDC method, based on the standard first-order IMEX Euler method, first for ordinary differential equations (ODE), and thereafter its application to the advectiondiffusion equation. Stationary domains are considered first, thereafter time-dependent domains.
In IMEX methods stiff terms are treated implicitly, in order to avoid prohibitively small timesteps [4]. The non-stiff terms are treated explicitly, as is often more computationally efficient.

IMEX spectral deferred correction as applied to an ODE
Here a shortened version of the methodology in [41] is presented. Consider the initial value problem for a first-order ODE on standard form where the solution U and a given initial value U a are in R, and F : R × R → R. It is assumed that F is smooth and can be decomposed as where F E is a non-stiff term and F I is a stiff term. The subscripts denote that the stiff part is treated implicitly and the non-stiff part explicitly. The interval [t a , t b ] is subdivided into p subintervals by choosing points t a = t 0 < t 1 < . . . < t p = t b , which are also referred to as instances (in time), and the length δt m = t m+1 − t m of a subinterval will be referred to as a substep. Given U a , initial (or provisional) approximate solutions U 0 m ≈ U (t m ) are computed at the instances {t m } p m=1 . To this end, introduce P = p + 1 Gauss-Lobatto points and use the standard first-order IMEX Euler method for these points. Then, the initial approximate solution is the solution to An increasingly accurate approximation U k+1 m p m=0 , for k = 0, 1, . . . K − 1, is obtained through a sequence of corrections. The superscript denotes the number of times it has been corrected; e.g. U k m is the kth corrected approximation of U (t m ). Each correction revises the current approximate solution and increases its formal order of accuracy.
A corrected approximation is computed as follows. Given U k m p m=0 we proceed as in [41] and obtain . To obtain a Kth-order method, it is sufficient to set P = K, i.e. p = K − 1, as explained below. Finally, we accept U K−1 K−1 as the approximation of U (t b ). In the right-hand side of (26) the expression I m+1 m U k m p m=0 appears, which is a pth-order numerical quadrature approximation of tm+1 tm F τ, U k (τ ) dτ. (27) Note that U k is a continuous function. This integral over [t m , t m+1 ] is computed by integrating the pth-degree polynomial that is interpolating U k over [t a , t b ], giving The coefficients q j m can be precomputed, reducing the quadrature to a simple matrix-vector multiplication. Each of the p different quadrature rules only require a rescaling of nodes and weights from the basic interval to [t a , t b ].
The error in computing the integral (27) with (28) is O δt P +1 , assuming F is smooth enough. Now, each application of the update equation picks up at least a single order of accuracy [7], thus to obtain a Kth-order method K − 1 corrections are sufficient. This implies that the number of Gauss-Lobatto points P should be set to be equal to K, as proposed in [41]. Thus, the total number of implicit systems, i.e. (25) and (26) To summarize. First we solve the given ODE (22) , given an initial value U 0 0 at t 0 , with a standard first-order IMEX method. Then the corrected terms U 1 are computed, initialized at U 1 0 = U 0 0 . This process is repeated K − 1 times. Finally we obtain the approximation

IMEX spectral deferred correction for the advection-diffusion equation on stationary domains
Consider the interval [t a , t b ] as a time-step [t n , t n+1 ] = [t n , t n + δt n ], which is subdivided as described above. Assume that the domain Ω is stationary. Writing the advection-diffusion equation (1)-(4) with the split into a stiff and a non-stiff term, as in (24), yields Inserted into (25), and collecting unknown terms on the left-hand side results in the modified Helmholtz equation (5), Recast as the modified Helmholtz equation (5) we have and (32) For both formulations of the modified Helmholtz equation we enforce the Neumann boundary conditions g, as evaluated at t m+1 .
Thus at each substep (m, k) both the potential ∇U k+1 m+1 and the Laplacian ∆U k+1 m+1 are computed. The former can be computed alongside u by (14) and (16), and the latter can be obtained directly from (5) once u is computed. These values have to be given as initial conditions as well at the start of each substep, as we need to initialize IMEX SDC with F E t 0 , U 0 0 and F I t 0 , U 0 0 for each time-step.
Remark 1 High order IMEX methods can also be used to compute provisional solutions in SDC, such as R-K methods. However, even though the solution in R-K methods is a linear combination of stage values, it is not necessarily true for its gradient and Laplacian. Some R-K methods have a property sometimes called first-same-as-last (FSAL), meaning that the final implicit stages are the same as the first ones in the next iteration, thus reducing the number of implicit problems to solve by one per time-step. The explicit stages are assumed to have closed forms, in which case FSAL is superfluous for the explicit part. In the method proposed in this paper, the explicit part contains −v · ∇U which is not explicit in a true sense; it is the gradient of the solution, which has to be solved for; or the differential operator has to be discretized, which is hard to do to high order. To the extent of our knowledge, only first-order IMEX Euler is FSAL also for the explicit term.

IMEX spectral deferred correction for the advection-diffusion equation on time-dependent domains
Here we discuss how to handle time-dependent and deformable domains when solving the advectiondiffusion equation (1)- (4). The problem we address is the following. Assume there is a timedependent domain Ω(t) superimposed on a discrete stationary grid X. When the domain Ω(t) changes over time there may be discrete points from X entering the domain where data from previous time-levels do not exist. These points are in the parts of the new domain that do not overlap the previous one, see figure 2. To solve the modified Helmholtz equation at this instance in time, all data need to be known at all discrete points inside the new domain. We take the approach of accommodating for the missing data by extrapolating existing data to these points using PUX, as described in section 4. First we present how we update the domain in time. We consider a Lagrangian specification of the flow field v. The position of a point x in a time-dependent domain Ω(t), or on its boundary Γ(t), changes in time as In our examples, the velocity field v will either be given explicitly, or obtained by numerical solution of the Stokes equations. In the latter case, the Stokes equations are solved with a boundary integral method, similar to the one presented in this paper, using the same discretization of the boundary. See the papers [42,5]. The update of x in (33) can be computed analytically in some simple cases, otherwise it is computed numerically, using an explicit R-K method of the same order as the IMEX SDC for the advection-diffusion equation. Clustering of the discretization nodes on the boundary is avoided by modifying the tangential velocity components [42]. Ω n Ω n+1 Figure 2: The domain Ω n (left circle) at time t n is advanced in time to Ω n+1 (right circle) at time t n+1 . The shaded blue area is the slice S n,n+1 where data in Ω n+1 is missing. Note that the size of the blue area is larger than for a typical time-step.
We consider an IMEX SDC applied to an interval [t n , t n+1 ], with the instances {t j } p j=0 , where t 0 = t n and t p = t n+1 , and δt m = t m+1 − t m for m = 0, . . . , p − 1. When addressing moving domains it is useful to see f from (32) as a sum over data for a given m and k. Thus the different f j correspond to different domains Ω j , with j = 0, . . . , p at t m for m = 0, . . . , p. Solving for U k+1 m+1 in (30) means solving the modified Helmholtz equation in Ω m+1 ∩ X, where only f m+1 is entirely known. The other terms f j are only partially known and need to be extrapolated to the points in Ω m+1 ∩ X where data does not exist, otherwise f will be incomplete. This is visualized in figure 3, which demonstrates the application of IMEX SDC of order three to an one-dimensional domain translated from the left to the right. Here we solve the governing equations inΩ 1 ∩ X, thus the data f 0 and f 2 need to be extended from Ω 0 ∩ X and Ω 2 ∩ X, see figure 3.
This motivates the definition of a slice S ta,t b . It consists of the points in X that are in the relative complement of Ω(t a ) in Then for each f j , with j = m + 1, no data exists in S tj ,tm+1 , we say that data is missing, and must be created artificially by extending it from Ω j ∩ X to S tj ,tm+1 . A slice is naturally thin, typically of the size of δx, thus the data is only extended a short distance.
We expect the procedure of extending data to the slices to impose a CFL-type condition of the form max The intuition is that the method becomes unstable if the slices are too large, i.e. the domain moves too much in a single time-step, relative the resolution of the underlying grid X.

Adaptivity
For an IMEX SDC it is by construction simple to obtain an approximation of the 2 error at t m after k + 1 corrections, namely r = U k+1 m − U k m 2 [11]. If the estimate r of the error is less Figure 3: A visualization of a three overlapping one-dimensional domains Ω 0 , Ω 1 and Ω 2 at different time-levels. Data is available at the black dots.
than some tolerance TOL , then U k+1 m is accepted as solution at time t n+1 . The time-step is then updated to δt n,NEW by δt n,NEW = δt n,OLD (0.9 TOL /r) with K being the order of the method and the value 0.9 is a safety factor. If the solution is accepted, then δt n+1 = δt n,NEW , otherwise the computations start over at t N with δt n = δt n,NEW . Thus even if the solution is accepted the step size is updated by the scheme (36), meaning that the time-step size may increase if appropriate.

Evaluating the particular solution
Recall the periodic formulation of the particular problem (11). It assumes the existence of an extension f e satisfying (12), that is f e = f in Ω and supp (f e ) ⊂ B, which is constructed with PUX as described in section 4. The approximate solution is expressed as a truncated Fourier series. If the data is represented on a uniform grid on a geometrically simple domain, then it allows for efficient use of FFTs to compute the Fourier coefficients. To this end, consider a box . These elements are referred to as nodes or points. The solution to a modified Helmholtz equation is computed at all points that fall inside Ω. It is straightforward to compute the particular solution (13) and its gradient (14) on X with standard FFT libraries. It involves two applications of a FFT, each of complexity O N 2 u log N 2 u . The order of convergence is q + 2 for the Fourier coefficients of f e ∈ C q (B) [17,53].
The gradient of the particular solution ∇u P is a correction in the boundary data for the homogeneous problem, see (10) and (19). It is evaluated on the boundary Γ with a nonuniform FFT (NUFFT) [21]. If N Γ is the number of evaluation points on the boundary, then computing the particular solution on these points by NUFFT is of complexity O N 2 u log N 2 u + N Γ , where N Γ is the number of discretization points on the boundary Γ, introduced in the next section.

Computing the homogeneous solution
Here we present the layout for solving the boundary integral equation (19) for the layer densities µ and how to compute the layer potentials (15) and (16) as a post-processing step.
Partition the boundary into N p segments, which we refer to as panels, such that For each panel there exists a parametrization and each panel is discretized in the parameter t with a N q -point Gauss-Legandre quadrature rule, using the canonical quadrature nodes {t G j } Nq j=1 and quadrature weights {ω G j } Nq j=1 . Then, a layer potential for some continuous kernel M over a single panel is approximated as where y ij = y(γ i (t G j )), µ ij is an approximation of µ(y ij ) and s ij = |γ i (t G j )|. Also, it will be useful to introducen ij =n x ij . If the geometry and integrand are well-resolved, then the approximation (39) is accurate to order 2N q . The kernels (15) and (16) are singular and nearly singular unless the target point x is well-separated from the source points {y j } Nq j=1 , resulting in a significant loss of accuracy. It can be regained by applying corrections obtained by product integration methods, as explained below.
We note that a parametrization of the boundary is only needed for initialization of the advection-diffusion equation solver. As the domain evolves in time it is the discretization that is updated; this applies to all methods we employ to solve the advection-diffusion equation in this paper. We use the approach in [46] to numerically update the discretization of the moving boundaries.
Discretizing (19) and collocating at the N Γ = N p N q nodes, known as the Nyström method, gives the linear system with D given by (20). As for the original boundary integral formulation, the Fredholm alternative states that there exists a unique solution. The corresponding system matrix is dense, but since its eigenvalues have only one accumulation point, which is away from the origin, the layer density µ can be solved for efficiently with GMRES [54]. The condition number and the dimensionality of the Krylov subspace are independent of the resolution of the boundary, assuming the density µ and the geometry are sufficiently resolved. However, the condition number does depend on the curvature [30,31].
After solving for the density µ, the solution u H and its gradient ∇u H can be computed by and with the kernels (17) and (18). Computing the sums in (40), (41) and (42) can be accelerated with spectral Ewald summation, see section 3.4. Then, the complexity is of order O(N Γ log N Γ ) for solving the boundary integral equation, and O(N 2 u log N 2 u ) to evaluate the solution on the uniform grid. Note that using standard Gauss-Legendre quadrature introduces errors due to the presence of near singularities when a target point is not well-separated from the source points. The numerical solution is corrected with kernel-split quadrature with product integration as explained below.

Kernel-split quadrature with product integration
Standard Gauss-Legendre quadrature methods are only accurate if sources and targets are wellseparated, otherwise the kernels (40), (41) and (42) are singular or nearly singular, which in turn degrades the accuracy in the approximate solution. To regain the lost accuracy we apply a kernelsplit product integration method for the modified Helmholtz equation, introduced in [19]. It is an extension of [26], which is a quadrature method based on decomposing the integrand into explicit singular terms multiplying smooth functions, and applying product integration; meaning that the smooth functions are approximated by polynomials and each term is integrated separately. This original method [26] is efficient and accurate for several PDEs, such as the Helmholtz equation and the Stokes equations [25,42]. However, as we explain later in this section, it may fail completely for the modified Helmholtz equation, hence the need for a modified version.
To apply the kernel-split quadrature in [19,26] the kernels (20), (17) and (18) must be split in order to identify the character of singularities. All three contain K 0 or K 1 which by standard decompositions [10, §10.31] can be written as where smooth remainder terms are collected in K S 0 and K S 1 . For on surface evaluation we have with which goes to zero in the limit as x goes to y. By (21) we have For S we have Here and we have no finite limit for S(y, x) corresponding to (47). Again the remaining terms are collected in S S , which is smooth. For the kernel ∇S the decomposition is where the symbol ∇ on the right-hand side is not the gradient, but part of the factors' names.
Here, ∇S S is a smooth function, ∇S C = 1, and Again, this has a log-type near singularity, but also a Cauchy-type near singularity.
In all three cases above the smooth terms are computed using regular Gauss-Legendre quadrature, but product integration is needed for all singular or nearly-singular evaluation in order to compute the layer potentials accurately. The method in its entirety is explained in [19].
We now demonstrate why product integration from [26] may be insufficient for the modified Helmholtz equation, and to that end only present a simple case with the log-kernel. Thus, consider a panel with end points −1 and 1, following the real line. The goal is to compute using the definition of the complex logarithm. Here, σ α is a smooth function that depends on the parameter α. Instead of applying standard quadrature methods, we expand the smooth function in a monomial basis. That is, and we have with the remainder R Nq = O(N q ). The involved integrals can be computed analytically, given by recursion formulas. Now let σ α (r) = I 0 (rα), as in (49), which grows asymptotically as exp (αr) / √ 2παr, where α is proportional to the inverse square root of the time-step. For large α, i.e. small time-steps, the function σ α may not be accurately represented by (53) with e.g. a 15th degree polynomial, or even a 31st degree polynomial, unless the maximum value of r is small. Consequently, product integration does not produce an accurate correction. A naive approach is global refinement of the boundary, however both K 0 and K 1 decay asymptotically as π/(αr)e −αr , i.e. they are very localized and only a small portion of the boundary needs to be refined. This would carry an unnecessary computational cost. Hence, the method in [26] needs to be modified.
In [19], new panels are created through recursive bisection such that they are of suitable length for product integration. It does not add degrees of freedom, but merely interpolates known quantities to shorther panels. The Cauchy-type singularity is also computed on the refined panels, as it increases the accuracy. Note that its multiplying function is a constant, thus product integration is not affected by increasing α.
The recursive bisection algorithm also addresses another problem: D L (y, x) log y − x and S L (y, x) log y − x grow with opposing signs to D S (y, x) and S S (y, x), both faster than K 0 and K 1 decrease. Thus they must cancel, potentially introducing cancellation errors when in finite precision.

Fast summation
Discretizing the integrals (19), (15) and (16) that arise from a boundary integral formulation yield the discrete sums (40), (41) and (42). Each can be computed directly at a cost O(N 2 ) for N sources and targets, which becomes costly for large N , e.g. when solving equation (40) with N = N p N q . In addition we need to compute several such sums in each time-step and inside each GMRES solve, where the number of iterations is independent of N . When solving the advection-diffusion equation with an IMEX SDC the O(N 2 ) cost can be decreased by fast summation methods, such as an FMM [6] or a spectral Ewald summation [1,40,38,2,45], to O(N ) or O (N log N ).
The references above are all for the kernels associated with more common PDEs, such as the Laplace equation and the Stokes equation. There is an FMM implementation for the modified Helmholtz equation described in [32], but to the best of our knowledge there are no Ewald summation formulas derived, upon which a spectral Ewald method can be built. In this paper we present a derivation of such formulas, along with truncation error estimates, and a spectral Ewald summation method for an O(N log N ) evaluation for the sums (40), (41) and (42), and their periodically extended counterparts. For the interest of space, the main results are presented in this paper, while supplementary material is found in [43].
We consider fast Ewald summation over an FMM since periodicity is simple to incorporate and considerable speedups are possible for uniformly distributed target points, such as the uniform grid, on which we evaluate the solution used in our computations. The basic concept of Ewald summation is as follows. A sum is decomposed into a sum in "real space" and a sum in "k space"; the terms in the former decays rapidly in R 2 and the Fourier transform of the latter converges rapidly in the frequency domain. Hence, the sums can be truncated at a cut-off radius and a maximum frequency, respectively. The evaluation of the Fourier part is then accelerated using FFTs, as will be explained below.

Ewald decomposition
We present the basic Ewald decomposition for u H and ∇u H in free-space for a target points x in Ω. The derivations are presented in appendix A. Initially we consider u H , the principle for ∇u H is analogous. Associate each pair of integers (i, j), i = 1, . . . , N p and j = 1, . . . , N q with some n = 1, . . . , N Γ , such that µ n = µ ij , s n = s ij etc., and write (41) as where f S,n = µ n s n ω G n .
By using the Ewald screening function the kernel S can be written as With this decomposition, an alternative formulation of the homogeneous solution u H (55) is with the closed form expressions S F (k, ξ) = 2π Here K 0 (·, ·) denotes the incomplete modified Bessel function of the second kind of zeroth order [23], defined as If a target point and a source point coincide, i.e. x = y n , then the self-interaction term is added to (59). Analogously, the decomposition ∇u H from (42) is with f ∇S,n = f S,n for n = 1, . . . , N Γ , and with K −1 (·, ·) from (62). The self-interaction term ∇S self is identically equal to zero. The sum in (40) can be treated analogously to ∇u H , since with f D,n = µ n s n ω G n . For the periodic case the sums over S R and ∇S R also include periodic images, and the inverse Fourier transforms are replaced by inverse Fourier series, see [43].

The spectral Ewald method
The spectral Ewald method is a fast summation method that makes the real space sum cost O(N ) and the k space sum cost O(N log N ) for N sources and targets [1,38,39]. The balance, in terms of computational effort, between the two is controlled by the splitting parameter ξ.
The summation in real space is accelerated as follows. Introduce a cut-off radius r c > 0, and include only contributions from source points within a distance r c of the target point x in the real space sum. That is with S R as in (60). The complexity can now be reduced to O(N ) by modifying r c with increasing N , such that the number of elements in {y n : y n − x < r c } remains constant. By first creating a linked list, the creation of {y n : y n − x < r c } also costs O(N ). The summation of ∇u H,R is treated analogously. We now study the k space sum u H,F . The goal is to compute it for any target point in the plane, discretizing the involved integral with the trapezoidal rule and accelerate the computations with FFTs. If α is small, which is not the case in this paper, then we need to take the approach discussed in appendix A.
The idea is now to spread the data at the source points y n to a uniform grid of N E points with spacing h on the box [−R/2, R/2] 2 , apply an FFT, scale the results with the Green's function and a window function, apply an inverse FFT and evaluate the integral in the k space sum at the target point x with the trapezoidal rule. The spreading means convolving the data f S with a window function w(x, ξ, η), where η > 0 is a scaling parameter. The Fourier transform of the window function is denoted w k , ξ, η := w k . We use for x = (x 1 , x 2 ). Here p denotes the number of grid points in the truncation, such that the domain of support is a square with sides of length ph. The parameter η balances the errors from truncating the window function and how well the Gaussians (70) are approximated on the grid. We use η = 0.95 2 πp/2. In [43] it is shown that if where then the k space sum approximation (124) can be written as whereH is the inverse Fourier transform of H . This integral is computed with the trapezoidal rule.
To summarize, the steps for computing the spectral Ewald approximation of the k space sum for u are • FFT : compute H using a 2D FFT, at a cost O(N E log N E ).
• Scaling: compute the tensor product H • IFFT : computeH(x) on the uniform grid using the 2D inverse FFT, at a cost of O(N E log N E ).
• Quadrature: The integral (73) needs to be computed for each target point x to obtain u H,F (x). This is done with the trapezoidal rule. Using the compactly supported window function (70), this is of cost O p 2 for each target point.
The steps for evaluating spectral Ewald approximation of the k space sum for ∇u are analogous, except now three inverse FFTs are required, since this is a vector quantity.
To further reduce the cost of the spreading and quadrature steps fast Gaussian gridding can be applied [21]. See [39] for implementation details.
For the sums (41) and (42) the target points x belong to the uniform grid X, and the final quadrature step can be avoided. It is achieved by constructing the uniform grid in the spreading step such that X is a subset of it. See [43] for how to modify the scheme accordingly.
To understand the computational complexity of the full scheme, we have to consider also the evaluation of the real space sum. We will comment on this in the parameter selection section 3.4.4 below.

Error estimates
The basic idea of Ewald decomposition is to split a slowly converging sum into two rapidly converging sums in their respective spaces. The resulting sums are then truncated based on carefully selected rules such that a certain tolerance is satisfied with optimal efficiency. These rules are presented here. For their derivations, see [43].
Truncation errors for the real space sums u H,R (68), and the corresponding for ∇u H,R , are introduced when they are limited to source points within a distance r c of the target point. Thus these errors are essentially the same in both the free-space and periodic setting. Let δu H,R and δ∇u H,R be the RMS errors. They can be approximated as and where Q f S = n f 2 S,n and Q f ∇S = n f 2 ∇S,n . How well the error estimates perform, for a range of values for ξ and r c , is demonstrated in figure 4 and in figure 5. They predict the errors well for error levels below 10 −2 .
For the k space sums u H,F and ∇u H,F truncation errors are introduced when wave numbers k that satisfy k ≥ k ∞ are excluded. The error estimates for the corresponding errors δu H,F and δ∇u H,F are different for the free-space and periodic setting. We present here only the former, the latter is given in [43]. The RMS errors can be approximated by and The truncation errors and estimates for u H are shown together in figure 4, for a test domain with 500 randomly distributed sources and targets, together with random point sources f S,n ∈ [0, 1]. In this example, α = 1 and L = 2π. The estimates hold for varying L and α. The estimates follow the error well for error levels below 10 −2 . For larger errors, the estimate is not as sharp. That is, however, far from the region of interest for an accurate computation. The corresponding plots for ∇u H are shown in figure 5, where we get similar results.

Parameter selection
The parameters of the Ewald summation impact the split of computational work between the real space sum and the k space sum and the level of the truncation errors. To keep the evaluation of the real space sum O(N ), for N sources/targets, the cut-off radius r c is set initially, such that the average number of neighbors of a target point is constant as N is increased. Using the truncation error estimate in (74) and (75) respectively, a suitable ξ is chosen given the desired truncation error level. Given ξ, an appropriate k ∞ = N E /2 is computed from the truncation error estimates for the k space sum (76), (77). This couples the grid size for the k space evaluations to N , the number of sources/targets, and yields a total computational complexity of the method of O (N log N ). The support of the window functions, p, is chosen large enough to keep the approximation errors low. Here, p = 24 suffices. The parameter η is chosen in order to balance the errors from truncating the Gaussians, and resolving them on the grid.

Partition of unity extension
PUX is a numerical method to construct high-order function extensions [17]. Discs, called partitions, are distributed along the boundary Γ; in each partition there is a Gaussian basis that interpolates f at the nodes in Ω that falls in the partition, then the interpolant is evaluated at the other nodes inside the partition, i.e. those not in Ω. Multiple partitions overlap and their respective extensions are weighted together with a high-order partition of unity to form a global extension. This has two purposes: to accommodate for missing data in the slices S n,m and to construct the extended function f e (12).

Interpolation with radial basis functions
We start out with a brief introduction of interpolation with RBFs and how to construct a local extension, a notion specified below. Recall that the computational domain is a box B = [− L 2 , L 2 ] 2 in R 2 that containsΩ. B has an underlying grid X consisting of N 2 u uniformly distributed points. Let ψ q be a compactly supported univariate symmetric positive function, such that ψ q (x) = ψ q ( x − c ψ ) for some point c ψ in X. We refer to this as a compactly supported radial basis function (RBF), centered at c ψ . The superscript q indicates the highest regularity subset C q 0 of C 0 that ψ q is a member of. The support of ψ q defines a partition ρ, i.e. ρ = supp(ψ q ) ⊂ R 2 , which is a disc with a radius denoted R. In the partition ρ there are N ϕ points {c ϕi } Nϕ i=1 drawn from the quasi-uniform Vogel-node distribution, defined as in local coordinates of the partition. This provides a near-optimal distribution of basis functions for interpolation [35]. Each point c ϕi is the center for a Gaussian ϕ i (x) = exp(−ε 2 x − c ϕi 2 ), where ε is a shape-parameter setting the width of the Gaussian. This forms a collection of Gaussians {ϕ i } Nϕ i=1 , distributed over the partition ρ, and it defines the interpolation basis in ρ. The standard form of RBF approximation I f of a function f in this basis is then The unknown coefficients {λ j } can be determined by collocating at the centers {c ϕi } of the RBFs and solving the resulting linear system The condition number of A cϕ is highly dependent on the shape parameter ε in the Gaussian basis: the smaller the shape parameter, to a limit, the smaller interpolation error, but the higher condition number; values of 10 18 are not uncommon for reasonable values of ε [34,13]. Also, the interpolation weights {λ j } are large in magnitude and are of varying sign. We can circumvent these drawbacks by proceeding as follows. Let X ρ denote the points {x ρ,i } in X within Euclidean distance R of c ψ . Furthermore, let f ρ be a vector corresponding to the values of f at X ρ . Then using (79) and Λ = A −1 cϕ f cϕ gives where A ρ,i,j = ϕ i (x ρ,j ) and A = A ρ A −1 cϕ . We now avoid using the interpolation weights, and are left with dealing with the ill-conditioning of the problem. It can be reduced significantly by using the algorithm RBF-QR [14]. The drawback is that RBF-QR is an O N 3 algorithm, and in PUX we consider several partitions; it is unfeasible to apply RBF-QR for every partition every time PUX is called. We return to how to circumvent this bottleneck later in this section.
We now describe how to create a local extension f e for the partition ρ. The setting is that f is only partially known in X ρ , and f cϕ is unknown and must be solved for. To do so, we introduce the following notation. Let the set X ρ be split into two disjoint sets: one set X I where function data f is known and X E where function data does not exist. Furthermore, let f e be a vector with unknown elements, corresponding to the unknown values of the function f in X E . We refer to them as the local extension of f . Let f ρ = [f f e ] where f is evaluated at X I , and suppose that both X I and X E are nonempty, and the interpolation matrix A is decomposed as then (80) takes the form We can obtain f cϕ by solving A I f cϕ = f in the least squares sense, as the basis locations and data locations are decoupled [35]. It is now straightforward to compute f e through The bottleneck of having to compute A with RBF-QR for every partition can be circumvented entirely, making it sufficient to call RBF-QR once as a precomputation step. It is based on the following observation. Assume that the number of Gaussian basis-functions is fixed, and that the center c ψ belongs to X, then A is independent of the position of c ψ ∈ X. It is also independent of the distribution of X ρ into X I and X E , as it merely affects the decomposition (81), i.e. the arrangement of the rows of the matrix A. Thus, RBF-QR is called once as a precomputation step to compute A, which is reused for all partitions for all time-steps.

Partition of unity
We now consider a collection of N ρ partitions {ρ i } Nρ i=1 , each associated with a compactly supported RBF ψ q i = ψ q ( · −c ψi ) where c ψi ∈ X for i = 1, . . . , N ρ . Thus all partitions are discs with radius R. Now, for each partition associate a weight function w i , defined as If ψ q i ∈ C q 0 , for i = 1, . . . , N ρ , then w i belongs to C q 0 . Also, if ψ q i is a positive function, then so is w i . The set of weights {w i } Nρ i=1 forms a partition of unity, meaning which is referred to in the literature as Shepard's method [49]. When a single partition is considered, we refer to the associated interpolation error, interpolant and so on as local, as opposed to global. Assume that there exists a local extension f e i for each partition ρ i . By weighting them together through a partition of unity we construct and we refer tof e as a global extension, or simply an extension, of f . This kind of extension has no imposed global regularity such thatf e ∈ C q (R 2 ); in most eventsf e will be discontinuous, as it is extended with zero outside of the partition. The global extensionf e is used for extending data into slices arising from time-dependent geometry, see section 3.1.3, as here, only data from within the partitions is used.

Function extension with compact support
We want an extension with a specified regularity as it is extended by zero outside its support, i.e. the global extension should have regularity q ≥ 0 in R 2 . To achieve this, the boundary of Ω is covered with overlapping partitions, i.e. the centers {c ψi } The weight functions (87) also form a partition of unity (85). The zero partitions are 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 = N ρ + 1, . . . , N ρ + N 0 ρ . Hence, as these zero valued functions are blended with the local extensions in the first layer of partitions, i.e. the extension partitions, the global extension will be forced to zero over the overlapping region. Therefore zero partitions should be distributed such that f e has a controlled transition to zero, Regularity ψ q (r) (1 − r) 5 + (8 + 40r + 48r 2 + 25r 3 + 5r 4 ) ψ 5 ∈ C 5 0 (1 − r) 6 + (6 + 36r + 82r 2 + 72r 3 + 30r 4 + 5r 5 ) Table 1: Wu functions ψ q ∈ C q 0 , with compact support in r ∈ (0, 1) [13]. Here (·) + = max(0, ·). The listed regularity excludes evaluation at the origin. see figure 6. Such a global extension f e of f is given by Here X Ω = X ∩ Ω and X E,i is X E for partition ρ i . Thus the only difference between how the extensions (86) and (88) are computed is how the weight functions are constructed. The global extension inherits the regularity of ψ q where it is extended to zero. However, there is no guarantee that f e , orf e , is of regularity q over the boundary of Ω. Inside we have f e = f , while outside f e is equal to an extension based on a weighted sum of local extensions. However, in practice we observe that the radial basis functions ψ q sets the regularity of the approximation.
As ψ q we use one of the compactly supported Wu functions, which are tabulated after their regularity q, see table 1 or [13]. 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. The partitions centers {c ψi } Nρ i=1 are set to be nodes on the regular grid that are the closest to be boundary, yet still in X Ω . 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. For instructions on how to set parameters to achieve high accuracy and a computational complexity linear in the number of grid points N u , as well as proofs for the regularity of the extension, see [17,16].
When the domains are time-dependent, both definition (86) and (88) are used for the global function extension. For each domain at each substep in [t n , t n+1 ] we construct an extension; first (86) is created for each domain which the modified Helmholtz equation is not currently solved on, i.e. create an extension without enforcing compact support with zero partitions. By doing so, missing data in the slices are created, which can be combined to create a right-hand side for the modified Helmholtz equation, see section 3.1.3. The right-hand side for the modified Helmholtz equation is in turn extended with (88) such that it has a prescribed regularity q in the plane.
Given two time stamps t j and t m+1 , using the notation from section (3. 1.3), available data f j in Ω j are extended into S tj ,tm+1 ∪B, meaning that for each time-step with SDC of order K, there are K − 1 + K(K − 1) 2 calls for PUX: there are K − 1 extensions for the initial approximations {U 0 m } K−1 m=1 . Then there are K − 1 applications of the correction equation, each with K values f j . The total cost is efficiently reduced by precomputing a single interpolation matrix A.

Implementation
To implement the software required to solve the advection-diffusion equation with the methods proposed in this paper we used MATLAB, except for the spectral Ewald method which is written in C++. There are also external libraries: for NUFFT we used the library from [20], based on [21]. The incomplete modified Bessel functions are not implemented in standard numerical libraries, but can be computed with the methods proposed in [24]; note that the precision is limited to ten digits. For the RBF-QR algorithm, used in PUX, we use the code from [33], as implemented in [14]. To label points as inside or outside Ω we use a spectral Ewald method for the Stokes potential [39], as well as compute the velocity field in numerical example in section 5.7 and section 5.8.
Let u * (x j ) and U * (t, x j ) denote approximate solutions to the modified Helmholtz equation (5), (6), and the advection-diffusion equation (1)-(4), respectively, at the point x j at the instance t, and let u and U denote the corresponding solutions or a computed reference solutions. Let e = u − u * or e = U − U * , then the 2 error and the ∞ error are defined as at time t, over the points {x j } NΩ j=1 which are the points in X that are in Ω(t).

The modified Helmholtz equation
To validate our proposed method for solving the modified Helmholtz equation (5), (6) with a = 0.3, b = 5, and c = −0.1045 + i5/439. If the Neumann boundary conditions and the forcing function are given by g = ∂u/∂n and f = α 2 u−∆u, respectively, then the exact solution is given by (90). For PUX we set the partition radius R = 0.15, which gives 77 extension partitions, and ε = 2; we keep this value for ε for all subsequent numerical experiments. An extension of f constructed with PUX, and a distribution of partitions along the boundary of the star-shaped domain, given by the parametrization (91), is shown in figure 6.
The boundary is discretized with N p = 200 panels, each with N q = 16 Gauss-Legendre nodes. For the first test we set α = 10 and solve with N u ranging from 40 to 1200. The results are shown in figure 7. All errors decay as O N −10 u , as expected [18,17]. Also, we see that the errors level out at 10 −10 and 10 −11 . This is due to the spectral Ewald method, where the involved algorithms for evaluating the incomplete modified Bessel functions are limited to about ten digits of accuracy [24].
The next test is to fix N u = 800 and study the errors' dependence on α. The errors for α ranging from 1 to 10 4 are shown in figure 7. For the relative ∞ error in u about one digit may be lost for α ∼ 10 3 , while for ∇u both errors clearly deteriorate, as several orders of magnitudes in accuracy are lost. Recall that α = O (Dδt) −1/2 , thus a high resolution in time results in more difficult problems in space. However, not until α ≈ 10 2 the ∞ error is larger than 10 −8 , for this particular problem, which corresponds to time-steps of magnitude 10 −4 ; the smallest value we consider for Dδt is of the magnitude 10 −4. 5 . From what we have observed this error does not limit the accuracy of U * when solving the advection-diffusion equation.

Convergence for rigid body in rotational flow
We now solve the advection-diffusion equation (1)-(4) on a rigid body in a rotational flow, in order to show third and fourth order convergence in time. Consider and a velocity field v given analytically by with w = 1, which is a positively oriented rotational field around the origin. Let (92) define the initial data U 0 at t = 0, the Neumann boundary condition g, and the source term J(t, x) = −2U (t, x), then (92) is an analytical solution. The domain Ω is given by the parametrization (91), with a = 0.05, b = 5, and c = 0 + i, see figure 8 for a visualization. The boundary is discretized with N p = 100 panels, each with N q = 16 Gauss-Legendre nodes. The periodic box B = [− L 2 , L 2 ], with L = 4.8, is discretized in each spatial dimension with N u = 500 uniformly distributed points; with these parameters the spatial error will not dominate the temporal error in the range of errors we investigate. The partition radius R in PUX is set to 0.15, giving 77 extension partitions.
We initialize the solver at t = 0, and measure the absolute 2 error at the terminal time t = 3, see figure 8 for a visualization of the solution at these instances in time. To show convergence we  successively reduce the magnitude of time-step δt. The results, using a third and fourth order SDC method, are shown in figure 9. We observe the expected third order convergence. We also observe the expected fourth order convergence, but only between two data points. Larger time-steps, without changing N u , violates a CFL-type condition, as we observe that the solution becomes unstable. At the other end the error levels out at 10 −10 as expected, due to limitations in accuracy in computing incomplete modified Bessel functions in the spectral Ewald method.

Conservation of mass
We study the conservative properties of our solver on a time-dependent geometry Ω(t) by computing the mass error e M (t) = M (t) − M (t 0 ) at time t, where Since Ω(t) cuts through the uniform grid, it is however difficult to numerically evaluate this integral in the general case, without introducing errors that obscure the errors in the mass conservation. Here, we do this for the special case where U * is the numerical solution to the advection-diffusion equation (1)- (4) with where k = 2π. Then the integral (94) is identically equal to 2π for U * = U for all t. The rotational field is given by (93) with w = 1, D = 10 −2 and the domain Ω(t) is the unit circle centered at c(t) = (− sin(t), cos(t)). The integral (94) is approximated by transforming to polar coordinates and then applying Gauss-Legendre quadrature in the radial direction and the trapezoidal rule in the angular direction, each with 150 nodes. We run the simulation from t 0 = 0 to the terminal time t = 3. The boundary is discretized with N p = 100 panels, each with N q = 16 Gauss-Legendre nodes.
We compute the mass error at t = 3 for a range of values for N δt = t/δt and N u , using SDC of order 3 and order 4. We also compute the max error e ∞ at the terminal time by comparing U *  with a reference solution using N u = N δt = 800 and SDC order 4. The results are shown in table 2. Clearly the mass error is at the level of the max error, or below. Furthermore, N u = N δt = 100 is clearly insufficient for SDC of order 3 to yield useful results.

Stability
Since we are treating the diffusion term implicitly and the advection and source term explicitly, we expect to have a first order CFL condition of the form (35). For time-dependent geometries, we extend data into slices S. Here, we want to investigate if this extension affects the CFL condition.
To do so, we solve the advection-diffusion equation on two different domains. Consider the unit circle centered at the origin with a star shaped domain inscribed in it. The star shaped domain is given by the parametrization (91) with the same parameters as in section 5.3, except that the radius is set to one and it is centered at the origin, see figure 9. The rotational velocity field (93) is imposed, thus the slices S are empty for the circle, but non-empty for the star shaped domain. We solve the advection-diffusion equation given in section 5.3 for a variety of values of w, δt, δx and D to study if the solution becomes unsteady for the star shaped domain, but not for the circle. We set L = 4.8 and run the simulation to time t = 3. The boundaries are discretized with 100 Gauss-Legendre panels, each with 16 nodes. The results are summarized in table 3, which shows that using PUX to accommodate for missing data, due to time-dependent geometry, does not impose a stricter CFL-type condition than the overall method. We leave further investigation of the CFL-type condition of the advection-diffusion equation solver to the future.

Adaptivity for rigid body in rotational flow
We now solve the advection-diffusion equation in the setting of section 5.3 with an adaptive time-stepper, as introduced in section 3.1.4. We run the code for the tolerances 10 −5 , 10 −6 , 10 −7 and 10 −8 for time-step orders K = 3 and K = 4; for the tolerance 10 −8 we set N u = 500, and N u = 400 for the other tolerances. The other parameters are set as in the previous section. We measure the error by comparing with the analytical solution (92) at each time-step, from t 0 = 0, to the terminal time t = 3. The absolute 2 errors over time are plotted in figure 10. We see that SDC for both orders K perform well, but K = 4 gives results closer to the set tolerance. With a Kth order method in total N mH = N δt K(K − 1) modified Helmholtz equations are solved for N δt time-steps, including the time-steps that were rejected due to not satisfying the set tolerance. A compilation of these values is presented in table 4. Clearly, K = 3 is preferable in terms of reducing N mH for the tested tolerances above 10 −8 . However, the errors are at some points about one order larger than the set tolerances, see figure 10, thus smaller time-steps would be required to satisfy the tolerance, and consequently increasing N mH . A safety factor in (36) less than 0.9 is recommended if using SDC of order K = 3 for less strict tolerances.    Table 4: Here we report data from solving the advection-diffusion equation for a rigid body in a rotational flow, see section 5.6, with an adaptive time-stepper of orders K, and different tolerances. The reported error is the maximum over the entire simulation. In the table N δt is the total number of time-steps, including the rejected time-steps, and N mH = N δt K(K − 1) is the total number of modified Helmholtz equations solved.
From studying convergence and adaptivity for a rigid body motion in section 5.3 and section 5.6 we conclude that, in terms of minimizing N mH , i.e. the total number of modified Helmholtz equations to solve, a fourth order SDC method is preferable over a third order SDC method only when considering tolerances stricter than 10 −7 .

Periodic channel
We now consider the advection-diffusion equation with homogeneous Neumann boundary conditions in a singly periodic domain, with an adaptive time-stepper as described in section 3.1.4. The domain Ω is a channel with sinusodial walls in a box B = [−L/2, L/2] where L = 4, that is periodic in the x 1 -direction. The initial data U 0 , with t 0 = 0, is which satisfies the imposed homogeneous Neumann boundary conditions, see figure 11 for a visualization of the initial data and the geometry. The source term J is set to be identically equal to zero for all t > 0, and D = 0.1. The walls are translated with a constant velocity of magnitude 4, from the right to the left. The velocity v is found by solving the Stokes equations in the domain with an imposed pressure gradient of (24, 0), i.e. from the left to the right, and with no-slip boundary conditions. The Stokes equations are solved using a boundary integral equation [5], using the same discretization of the walls as for the modified Helmholtz equation. Using the periodicity of the problem, the velocity v(x) can be computed at every point in X ∩ Ω(t) for any value for t using post-processing. See figure 11 for a visualization of the flow field. We set N u = 400, and N p = 100 panels on each wall with N q = 16, and partition radius R = 0.3. A reference solution is computed with N u = 400 and a tolerance 10 −10 for the estimated absolute 2 error, with K = 4. The absolute 2 error is measured at t = 1 by comparing the numerical solution with the reference solution. The results are presented in table 5. Two snapshots of the numerical solution at t = 0.5 and t = 1 are shown in figure 12.
For K = 4 and tolerance 10 −8 the absolute 2 error is slightly above the set error tolerance; for the other tested less strict tolerances the 2 error is well below the set tolerance. For K = 3  we get similar results for the tolerance 10 −8 , but by solving 1230 more modified Helmholtz equations. As for the numerical example in section 5.6, the error is one order larger than the set tolerance for tolerances 10 −5 and 10 −6 . The results are consistent, as we observe again that K = 3 is preferable, in terms of efficiency, over K = 4 for less strict tolerances.

Deforming drop
We now consider the advection-diffusion equation on a time-dependent domain Ω which deforms under the effect of a given velocity field, resembling a deforming drop. This is done for two different initial values and velocity fields: in the first setting we show convergence in time, and in the second we demonstrate robustness. Consider the initial data (96) at t = 0 in a circle of radius 1.5 centered at the origin, where ∂U 0 /∂n = 0 on the boundary. As in the previous example we impose homogeneous Neumann boundary conditions and set the source term J identically equal to zero. The diffusion coefficient D is set to be equal to 0.1.
The time-dependent velocity field is given by v(x) = (cos(πt) cos(x 2 ) sin(x 1 ), − cos(πt) sin(x 2 ) cos(x 1 )), which is an extensional flow, meaning the circle at time t = 0 will become elongated and then go back to its initial state, see figure 13. After t = 0 we have no parametrization to describe the boundary, thus we move the discretization of it according to the velocity field. A reference solution is computed at time t = 1 on a box with L = 5, N u = 500, and δt = 3.1250×10 −4 . The boundary is discretized with N p = 100 panels, each with N q = 16 Gauss-Legendre nodes. The partition radius in PUX is set to 0.4, and the partitions are redistributed along the boundary as it changes over time.
The time-step δt is successively increased with a factor of two and we solve the advectiondiffusion equation with SDC of order 3 and 4. The absolute 2 error is computed at t = 1 by comparing with the reference solution. The results are shown in figure 13, where we clearly achieve third order convergence. Fourth order convergence is almost obtained; some other error or errors are interfering with the temporal convergence. For larger time-steps a CFL-type condition is violated, thus the convergence line cannot be continued. For smaller time-steps the error line would level out, again due to limitations in computing the incomplete modified Helmholtz equations. As seen before the third order SDC is more efficient for less strict tolerances. We now study the robustness of the method simulating a drop with more complicated initial data and velocity field. Thus let where k = 2π/1.5, c = (0, 0.5), in a circle with radius 1.5 centered at c. The imposed velocity field is a superposition of the a rotational field (93) with w = 1 and the extensional field (97). We run the code for the time interval [0, 5] for homogeneous Neumann boundary conditions and a zero valued source term J, with N u = 300, with the L and boundary discretization and partition radius as above. The time-step is 0.01.
The results are shown in figure 14. Clearly the method is robust and can handle deformations into entirely new shapes. As expected the solution goes towards an equilibrium state due to diffusion.

Conclusions
We have presented a boundary integral equation-based method for solving the isotropic advectiondiffusion equation in two dimensions on time-dependent geometry. First, the PDE is discretized in time by applying an IMEX SDC method, then an approximate solution in time is advanced by solving a sequence of modified Helmholtz equations. By using PUX [17] the necessary quantities are smoothly extended to accommodate for missing data due to time-dependent geometry. The data is represented on a stationary underlying grid, while a discretization of the boundary is updated in time. We obtain a tenth-order accurate method in the resolution of the underlying stationary uniform grid when solving the modified Helmholtz equation. To reduce the computational complexity of computing the involved sums in the boundary integral method we present Figure 14: Snapshots of simulation based on numerical solution of the advection-diffusion equation with initial data (98) at different times. Notice the different ranges of the color bars. From left to right, top to bottom: t = 0, 0.5, 1.5, 2.5, 3.5, 4.5. a spectral Ewald method, which is of complexity O(N log N ) for N targets and sources. The time-stepping method is adaptive, providing errors down to 10 −8 , and we show up to fourth order convergence in time for time-dependent deformable geometry. Moreover, the single interpolation matrix used in PUX is precomputed once and reused for all time-steps. In the precomputation step the condition number of the extension problem is significantly reduced.
This paper is an intermediate step in developing a complete solver for the simulation of drop dynamics in Stokes flow, with surfactants dissolved in the bulk fluid and on the drops' surfaces. This would mean full coupling between the advection-diffusion equation and Stokes equations.
One of the challenges that remains is to make PUX adaptive, as high resolution will be required in the gap between two closely interacting drops. Since PUX currently requires a uniform grid for efficiency, high grid resolution is imposed globally. Another future direction is to investigate alternatives to IMEX SDC methods, in order to reduce the number of modified Helmholtz equations to solve per time-step. Future work also includes closer investigation of the CFL-type conditions that are imposed by using elliptic marching with SDC.

A.2 Derivation of Ewald decomposition for ∇u H
We present the derivation of the Ewald decomposition for u H in the free-space case, seeking ∇S R and ∇S F in (64), and ∇S self . The derivation for the periodic setting is analogous and can be found in [43].
A. 3 Computing u H,F in the spectral Ewald summation for small α If the parameter α in the modified Helmholtz equation (5), (6) is small, then computing u H,F (x) in (59) with the trapezoidal rule, accelerated with FFTs, will not be accurate. This is due to S F containing the factor 1/(α 2 + k ) and therefore is nearly singular at k = 0. To rectify this, we proceed as in [2], based on the techniques in [55], and cut off interactions in physical space beyond the domain of interest, i.e. B. Introduce where R is larger than the largest point-to-point distance in the domain B, and rect(x) = 1, x ≤ 1/2, 0, x > 1/2.
(125) Note that when α 2 is large enough there is no need to modify S F and ∇S F , see figure 6 in [43]. This is the case when solving the advection-diffusion equation in this paper.