Fast Ewald summation for Green's functions of Stokes flow in a half-space

Recently, Gimbutas et al derived an elegant representation for the Green's functions of Stokes flow in a half-space. We present a fast summation method for sums involving these half-space Green's functions (stokeslets, stresslets and rotlets) that consolidates and builds on the work by Klinteberg et al for the corresponding free-space Green's functions. The fast method is based on two main ingredients: The Ewald decomposition and subsequent use of FFTs. The Ewald decomposition recasts the sum into a sum of two exponentially decaying series: one in real-space (short-range interactions) and one in Fourier-space (long-range interactions) with the convergence of each series controlled by a common parameter. The evaluation of short-range interactions is accelerated by restricting computations to neighbours within a specified distance, while the use of FFTs accelerates the computations in Fourier-space thus accelerating the overall sum. We demonstrate that while the method incurs extra costs for the half-space in comparison to the free-space evaluation, greater computational savings is also achieved when compared to their respective direct sums.


Introduction
The determination of the motion of particles in bounded or unbounded flows is a central problem in microhydrodynamics. For a large class of industrial processes like particle filtration, sedimentation or aggregation, and deposition of pulp fibres in paper manufacturing, the fluid inertia is negligible and the governing equations are well-approximated by the Stokes equations [3].
The system of Stokes equations is linear and can be reformulated as an integral equation. In a boundary integral method, once the integrals are discretized, discrete sums with fundamental solutions of Stokes flow remain. Typically, any exterior solid boundaries or interfaces between different fluids are discretized. Periodicity in one or more directions is however usually built into the definition of the fundamental solution, leading to an (infinite) summation of periodic images in the discrete sums. Thus one can speak of 1-periodic, 2-periodic, 3-periodic or free-space problems -indicating the periodicity built into the evaluation of the Stokes potentials. The finite free-space sums for evaluating Stokes potentials are of the form (or some variant of) (4), i.e., each summand is a convolution of the Green's function with a source term.
For a simple geometry like a flat plane in unbounded space, it is also possible to avoid explicit discretization of this plane, and instead modify the evaluation of the Stokes potentials to achieve a no-slip condition. This requires the introduction of sources at image locations reflected in the wall, as well as correction terms. Such an explicit representation was first derived by Blake [4] and later in a more elegant form by Gimbutas et al [1]. The fast evaluation of these sums will be the focus of this paper, and we start by giving some background to the problem.
If the Green's function is the harmonic kernel and the source term a scalar, the sum corresponds to the formula for the Coulomb potential of a system of point charges. The interest in a fast evaluation of such sums actually stemmed from this particular problem, with Ewald's investigation [5] of the 3-periodic case in 1921 now known as the Ewald summation technique. Hasimoto [6] then considered the 3-periodic sum of stokeslets. In these decompositions, some specific choices are made to turn a 2 The free-space problem The three fundamental solutions for Stokes flow are the so-called stokeslet, stresslet and rotlet singularities; they are tensors that have an explicit representation as follows: with the short hand notation r =: r and x, y ∈ R 3 . Here ⊗ is the standard tensor product for vectors in R 3 , I is the second order identity tensor in R 3 and E is the alternating third-order tensor whose representation in the natural basis of R 3 coincides with that of Levi-Cevita's symbol.
In the Einstein summation convention or index notation, the direct representations above reduce to 8πT ijk = − 6 r i r j r k r 5 , (2b) where δ ij , ijk are the Kronecker delta and Levi-Cevita's symbol respectively. Wherever possible we shall use the direct notation for economy and provide the index notation as explanation. Given that there are N point forces of intensity 8πf (m) at locations y (m) , m = 1 . . . N , the stokeslet induces a velocity field at x ∈ R 3 . Similarly, if there are N "double forces" of intensity 8πg (m) and orientation q (m) , the stresslet and rotlet tensors also induce a velocity as follows: where S (m) = S(r (m) ), and T (m) and W (m) are defined similar to S (m) with the short hand notation x − y (m) =: r (m) . The formula for the rotlet can also be written in terms of the vector cross product by replacing the skew-symmetric tensor with the axial vector as has been done in [2], but we do not do so in order to stay consistent with the notation in Gimbutas et al [1]. The case where the evaluation point x = y (p) , p = 1, 2, . . . N, p = m will be the one considered here on as it is of most interest and occurs in boundary integral methods and potential methods for Stokes equations. In index notation, the velocity field induced in each case is written as The result of summation over all N terms above in (4) gives the free-space velocity corresponding to a collection of stokeslets, stresslets or rotlets. Note that the operation in each case is a convolution 1 between the kernel and the source term. In [2], a fast SE method is proposed for this evaluation.

The half-space problem
A natural extension of this question would be to ask if there exists an explicit representation for the velocity field due to discrete singularities in the half-space, and secondly, if there are fast methods to evaluate it. This paper deals with the latter question. Indeed, there are closed form expressions for the half-space, those given by Blake [4] and Gimbutas et al [1] which use the method of images to derive an appropriate expression. While the older representation of Blake requires evaluation of multiple harmonic and dipole fields, the recent work by Gimbutas et al provides a very elegant representation using correction terms expressed in terms of a single harmonic potential. Their key idea was to invoke the Papkovitch-Neuber [20,21] representation formula for constructing a divergence-free velocity field that satisfies the Stokes equation using a harmonic potential. The free-space formula itself does not satisfy the boundary conditions at the wall, but the combination with its image ensures the tangential component satisfies the boundary condition. The Papkovitch-Neuber correction term is added to adjust the normal component of the velocity at the wall. The result is the formula (5) which is divergence-free, satisfies Stokes equations, and the no-slip boundary conditions at the wall. In this paper, we show how the SE method illustrated in [2] can be applied to this representation formula for the half-space (5) to yield a fast method of evaluation.
y (1) y (1) y (1) y (2) y (3) x 3 x 3 x 2 First we set out our notation in order to restate the formulae for the half-space. A system of Cartesian coordinates is set up such that the x 3 = 0 plane coincides with the wall.
Then the representation for the half-space in each case involves a harmonic function associated with each mirror image as follows: The harmonic functions φ S , φ T and φ W presented in [1] are written in terms of potentials due to point sources, dipoles and quadrupoles. Here we instead express them as gradients of the harmonic potential. To explain, suppose the harmonic potential G of a unit charge 4πG := 1 r , we can write the potential of a dipole with orientation a as G D := −∇G · a and that of a quadrupole with intensity b and orientation a as G Q := ∇∇Ga · b. Thus, using these relations the expressions for φ S , φ T and φ W recorded in [1] are restated as where G (m) = G(r (m) ).
Note that the formulae (5) are not translation-invariant, and it is essential that the origin of the coordinate system be located on the wall or boundary of the half-space.

The Ewald decomposition
Here we quickly summarize the motivation behind the Ewald decomposition. The idea is to introduce a scalar (Ewald) parameter ξ and split the fundamental solution or Green's function (for the stokeslet, say) into where S F is the Fourier transform of S F and F −1 indicates the inverse Fourier transform operator (IFT). The formula for u(x) then becomes The first term represents local interactions and it will be seen below that it decays exponentially fast, and can hence be truncated. Thus it is evaluated directly by summation over all sources located within a chosen distance r c of the target x. The second term is the IFT integral and if the integrand were smooth and compactly supported, the integral could be approximated to spectral accuracy with a trapezoidal rule in each coordinate direction, allowing for the use of FFTs for the evaluation. However, we will find that all the kernels and correction functions relevant to this present work will have a factor (like B(k) in the expression for S F (k, ξ) recorded below) which makes them singular at k = 0. The method to circumvent this singularity in our quadrature is one of the key points and will be discussed in detail subsequently.
If the target location x coincides with a source location, then the so-called self-interaction term must also be accounted for and removed. This is evaluated as the limit of S R − S when r (m) = x − y (m) → 0.
For purpose of illustration, we explicitly show the Ewald decomposition of the stokeslet derived in [6].
The corresponding expressions for the stresslet and rotlet are tabulated for reference in the appendix in equations (14) and (15).
The computational complexity of the direct sum in (4) is O(N 2 ) for x = y (p) , p = 1, 2, . . . N, p = m, and the Ewald sum by itself does not reduce the complexity. The SE method keeps the real-space sum at O(N ) by a specific choice/scaling of parameters when scaling up the system, and this combined with the use of FFTs reduces the over-all complexity (see [2] for details) to O(N log N ).
One would like to use the SE method for (5) also to accelerate the process. While at first glance, the expressions involved in the sum may appear to be too cumbersome and complicated, a careful perusal will reveal that there is a structure to the formula; the summand consists of a linear combination of a stokeslet/stresslet/rotlet and its image about the x 3 = 0 plane along with two other terms involving the respective harmonic functions, which could be viewed as corrections to satisfy the boundary conditions. The Ewald decomposition of the first two terms on the right hand side of (5) follows directly from [2], so it only remains to derive an appropriate Ewald decomposition for the two correction terms involving the harmonic function and its gradient. We recognize that the Ewald decomposition of G, ∇G, ∇∇G and ∇∇∇G is required to complete this task.
Turning to it, a lengthy (but straight-forward) calculation finds the real space components (denoted by the subscript R) to be where, r i is the i th component of the vector joining the source-location to the target-location, and f (r) := erf(ξr) r for convenience. The components evaluated in Fourier-space through their Fourier Transforms and denoted by the subscript F are where k = (k 1 , k 2 , k 3 ) T , k := k ,k i := k i k and H(k) := 4π k 2 . The reason for this notation will become clear later.
These correction terms are centered at the image locations but the evaluation point x is never an image location hence there is no need to account for the so-called self-interaction term.

The half-space Spectral Ewald method
The half-space formula described in Section 3 essentially transforms the problem into a new free-space problem, albeit with extra terms. Hence the computational framework used is the same as described in [2], following the recent idea introduced by Vico et al . [15] to solve free space problems by FFTs on uniform grids. Thus we keep our explanation of the method very general and brief, directing readers to [2] for complete details while we highlight only the differences from it and finer points of note. A simple but useful observation is that for computations, the first two terms (kernel and its image) can be combined into a single term by absorbing the negative sign into the image source vectors, so that one can think of a kernel with 2N sources.

The real space component
The real space components for the stresslet and rotlet derived in [2,12] are reproduced in Appendix 10, while that of the respective correction functions φ S , φ T , φ W and their gradients can be written down directly by utilizing (7) in (6).
A cell-list of nearest neighbours within the cut-off radius r c is prepared for each target, the difference with [2] being that all target locations were also source locations in a free-space evaluation, but for the case of a half-space, target locations are on one side of the wall and also serve as sources, but their reflections act as sources only and are located on the other side of the wall. Other than that, the procedure to evaluate the local interactions is as before. The choice of the cut-off radius r c is made from the desired error-level using the truncation error estimates discussed in Section 6.

The Fourier-space component
The Fourier-space component is evaluated through the integral which always has the form where the quantities A(k, ξ),K(k) depend on the choice of kernel (see Appendix 10,K is eitherB(k) orĤ(k)), while the quantity c (m) is a scalar, vector or tensor that depends on the m th source term and location. For ease of explanation we write out in full the expression that emerges for the stokeslet from (5), Substituting from (6) and (8), it expands to where we have rewritten the first integral as a sum over 2N terms by setting While the equation above has been written out in full for the stokeslet, the method that will be discussed carries over without modification for the stresslet and rotlet as well since the expanded formula in that case too has the same structure and number of terms. The explicit presentation of the five integrals that need evaluation on the right hand side of (11) serves a dual purpose: (1) It shows that the integrals arising from the correction terms have the same form as that of the integral arising from the stokeslet in free-space. For a single integral of the type considered here, Klinteberg et al [2] in Section 5.1 illustrate and justify the sequence of operations followed. In the case of the half-space however, we have five integrals of that type, and we shall now explain how to combine them together in the evaluation while following the same procedure.
(2) It also makes clear that all the integrands have the factor B(k) or H(k) which make it singular at k = 0.
To circumvent this issue, we will continue to follow [2] and [15] and introduce modified Green's functions where the Fourier transforms of these functions have no singularity at k = 0. The truncated Green's functions are denoted by the superscript R that stands for the radius of the support in real space, and their Fourier transforms are given by Using these truncated Green's functions will yield exactly the same result for the harmonic/biharmonic equation as the original ones, as long as the right hand side has compact support within the solution domain, and R is chosen sufficiently large. Specifically, if the solution domain is such that the largest point to point distance is R max , then we need R ≥ R max .
We assume that the physical domain is a cube of size L that encloses all the sources, the computational domain is a cuboid of dimensions L × L × 2L that encloses both the sources and their images.
In the Ewald decomposition, sources are convolved with Gaussians or modified Gaussians to form the right hand side that defines the Fourier space problem. Hence, the assumption above of a compactly supported right hand side is violated. In the actual discretization, we interpolate point sources to the grid using a window function. In the SE method, this window function is a (suitably scaled, hence not the same) Gaussian, as will be explicitly defined in Algorithm 1. The domain length L will be extended toL to accommodate the support around the source locations. With the parameter choices that we will soon detail, this extension is sufficient and further extension will not reduce the total error. A more detailed discussion of this issue can be found in section 5.3 of [2] . The resulting computational domain is discretized by an equi-spaced grid with spacing h containingM ×M × 2M = 2M 3 points. Note that such a grid induces, in its k-space counterpart, a spacing ∆k = 2π/L. Fundamentally, we are performing an aperiodic convolution in Fourier space, and hence need a twice over-sampled representation of H R (k), i.e., a representation with k-space resolution of ∆k/2 . It has been shown earlier by Vico et al [15] that the terms H R (k), B R (k) can be evaluated knowing only the value of R, which itself is determined by the size or extent of the domain. Thus for computational efficiency, we precompute H R (k), and B R (k) for the stokeslet or stresslet on a grid with 16M 3 points. We cannot compute this directly by starting from values on the physical grid since it is only the Fourier transform that is known analytically. This computation is thus carried out as follows: 1. Evaluate H R (k) and B R (k) on a grid of spacing ∆k/s f (or 2(s fM ) 3 points) where the truncation radius for the domain R = √ 6L and the oversampling factor s f ≥ 1 + √ 6 is chosen as small as possible such that s fM is an even integer.
2. Compute the 3D-IFFT and truncate to get H R on a grid of 16M 3 points.
3. Compute the 3D FFT now to get back H R (k) on a grid with spacing ∆k/2, that is a twice oversampled representation.
This set of values is now used in the algorithm. Note that this is different from simply sampling H R (k) on a 16M 3 grid. The reason becomes clear if we consider the formulae in section 4.3 in [2]. It is apparent that we need to truncate the Green's function values centered at a particular point in the physical grid. Thus, to perform the truncation, we start with sampling values of H R (k), perform an IFFT to obtain values in physical space, and then perform an FFT after truncation.
Algorithm 1 SE method for Fourier-space calculation of stokeslet for half-space 3 f (m) I play the role of c (m) in the general formula (9). 6: Compute the FFTs, i.e., C 1 (k), C 2 (k), C 3 (k), zero-padded by a factor of 2.
6L oversampled by factor of 2 (see remarks on precomputation). This corresponds to K(k) in (9). 8: Compute the product A(k, ξ) K(k)e −(1−η) k 2 4ξ 2 C(k) with reference to the general representation (9) appropriate for each of the 5 integrals, and where C(k) hence is one of C i (k), i = 1, 2, 3. Call them t 1 , t 2 , t 3 , t 4 , t 5 The first three are vectors and the last two are scalars. 9: Compute the IFFT of t 1 + (0, 0, t 4 + t 5 ) T and t 2 + t 3 , truncate to get T 1 , T 2 respectively oñ M ×M × 2M grid. 10: By using trapezoidal rule and truncated Gaussians at the target point x = y (m) , compute In comparison to the evaluation for free-space in [2], the number of FFTs and IFFTs increases and this information is summarized in Table 1. The increase in the number of FFTs is consistent with the half-space formula (5); the correction terms require one vector and one scalar FFT (3 + 1) for the stokeslet, one vector FFT (3) for the rotlet, and one tensor and one vector FFT (9 + 3) for the stresslet. The increase in the number of FFTs required is thus not surprising, but the increase in the number of IFFTs is not so obvious, for after calculating the FFTs and performing the convolution as a product in Fourier space, it might seem that we could combine them all and perform a single IFFT and integration step. However, that is not possible due to the presence of the coordinate x 3 multiplying the gradient. This term needs to be treated separately and hence it requires an extra IFFT and integration step. One might well ask if the Ewald summation technique is still worthwhile with this increase in the complexity of the problem in Fourier space; fortunately, the answer is affirmative, as demonstrated in the next section.

Truncation Errors
The errors in the real-space calculation are caused by truncation. However, the errors in the SE method for the Fourier-space part are not due to truncation alone; there are approximation errors due to the quadrature rule for integration and the discretization and truncation of the Gaussians in the spreading and quadrature steps. Given P , the number of points across a truncated Gaussian, the parameter η = η(P ) can be chosen to balance discretization and truncation errors, which leaves P as the only free parameter. Approximation errors decay exponentially with P and by using a sufficiently large P , the approximation errors may be considered negligible so that the measured error is due to truncation errors only [22], as the Ewald parameter ξ does not introduce any errors.
The truncation errors appear in the real-space part because we consider only local interactions within the radius r c while in the Fourier-space part, the integral in (9) over all of R 3 is truncated in practice to consider some large but finite wave number with magnitude k ∞ . The value of k ∞ is related to the grid and computational domain by the relation k ∞ = π/h = πM /L. The exact real-space contribution is obtained by letting r c → ∞ in the real-space sum, and in combination with the naive direct sum, it also yields the exact Fourier-space contribution 2 . From these, the computed truncation errors are found. In [2], the authors report truncation error estimates based on the methodology introduced by Kolafa & Perram [23] that agree closely with the computed errors for both real-space and Fourier-space evaluations. These are statistical error estimates for the root mean square (RMS) truncation error, defined as Implicit in these estimates are the assumptions that the sources are randomly distributed, and that the error measure has a Gaussian distribution.
Since the derived estimates are statistical in nature, the methodology of Kolafa & Perram will not be able to account for possible cancellations due to symmetry in the half-space formula. On the contrary, this approach will yield an estimate that is the sum of two contributions-1. Truncation errors due to the stokeslet, stresslet or rotlet kernels (as in the Tables)

Truncation errors due to the correction terms
Such an estimate is likely to be a conservative upper-bound for the computed errors. Therefore, we examine the truncation error expressions to ascertain if one can neglect the error contribution made by the correction terms.
From Tables 2 & 3, it is evident that the kernel estimates have an exponential decay term (e −ξ 2 r 2 c or e −k 2 ∞ /4ξ 2 ) multiplied by some power of r c or k ∞ . In case of the half-space, the correction terms are all harmonic functions and/or their derivatives, and their truncation error contribution (see [23]) decays much faster than the stokeslet or stresslet kernels and at the same rate as the rotlet kernel. Thus, the existing truncation error estimates for free-space are a good starting point for the half-space as well. Of course, in evaluating the estimate, the sum is over both sources and images, and this modifies the quantity Q and the RMS error.
In Figures 2 and 3, the truncation error estimate is compared to the computed error for both real-space and Fourier-space parts. This is done by calculating the relative RMS error which is the ratio of the RMS error and the RMS value of the velocity considered at all targets. For reference, the corresponding curves for the free-space problem are also plotted. For both the stokeslet and stresslet in Figures 2 and 3, the computed errors are lesser than the half-space estimate, and the freespace computed error and estimate. The computed error for the half-space here is lower due to the cancellations induced by symmetry. For the rotlet, while all these curves are much closer to each other, closer examination reveals that the previous trend is no longer upheld. This is because the correction terms that have been neglected have the same order as the kernel. However, the overall agreement of the computed error with the existing estimates justifies the decision to neglect the contribution of the correction terms to the truncation error estimate. Table 2: Fourier-space truncation errors for the stokeslet, stresslet, and rotlet [12] for free-space. The quantity Q is defined for each kernel as in the second row and · F denotes the Frobenius norm for second-order tensors. Table 3: Real-space truncation errors for the stokeslet [13], stresslet, and rotlet [12] for free-space.
The quantity Q is defined in the same way as for the Fourier component in the previous table.  Table 3. Red coloured dots and lines are associated with the half-space while blue colour is associated with the free-space evaluation.The system is N = 2000 randomly distributed point sources in a cube with sides L = 3, with ξ = 4.67, and r c ∈ [0, L/2].

Numerical Results
We consider N random point sources drawn from a uniform distribution from a box of dimension L × L × L. The sum (4) is evaluated with stokeslets, stresslets and rotlets, at the same N target locations. All components of the source strengths and source orientations are random numbers from a uniform distribution on [−1, 1]. All computationally intensive routines are written in C and are called from Matlab using MEX interfaces. The results are obtained on a laptop workstation with an Intel Core i7-6500U Processor (2.50 GHz) and 16 GB of memory, running two cores unless stated so. Actual errors are measured by comparing the result with evaluating the sum by direct summation. For a given system of N charges in a box of edge length L, the required parameters for our halfspace Ewald method are the Ewald parameter ξ, the real space truncation radius r c , the number of grid points M covering the computational domain, and the Gaussian support width P . Other parameters like δ L ,L,M are then set automatically from these (see Algorithm 1 and [2] for details). For a large-scale numerical computation, ξ, r c , M and P must be set optimally. For a given value of ξ and absolute error tolerance for the free-space evaluation, close-to optimal values for M and r c were computed in [2] using the truncation error estimates in Tables 2 and 3. We use the very same optimal values as starting points for the half-space since the same truncation error estimates hold good for the half-space too as shown. Then we perturb ξ to achieve the smallest runtime while keeping ξr c and k∞ ξ constant. Note that k ∞ = πM L . The optimal set of values found is used for larger systems by scaling k ∞ such that L k∞ is constant. The results for free-space in [2] convincingly demonstrated the need for the Ewald decomposition by exhibiting the speed-up gained over the naive direct sum. The present work aims to make a similar case for the half-space. Before that however, it would be interesting to compare the computational expense of evaluating the direct sums themselves for the free-space and the half-space formulae. This has been done in Figure 4 and it is seen that on average, the half-space sum is between 3.3-3.7 times more expensive (with rotlet the least, and stresslet the most) than the free-space sum and this factor stays constant even as the number of sources increases. The obtained range 3.3-3.7 is reasonable since the direct sum for the half-space involves 4N terms (kernel and correction terms), in contrast to that for the free-space which involves only N terms. The additional expense due to the correction terms is smallest for the rotlet, and most for the stresslet, and this is not surprising seeing the formulae (6). We now tackle the question of whether it is worthwhile to consider the Spectral Ewald method for  the half-space formula despite the substantial increase in the number of FFTs and IFFTs that need to be performed for the Fourier-space component. In Figure 5a the computing time for evaluation of the sums is plotted versus N , for all three kernels and for both the Half Space Ewald (HSE) method and direct summation. The relative RMS error for HSE is kept below 0.5 × 10 −8 . As we vary N , we maintain a constant density N/L 3 = 2500 by changing the size of the box. The system is thus scaled while other parameters of the method like ξ, r c , P and the grid resolution L/M are kept constant. As both N and L increase, so does the grid size. For all kernels, we set P = 16, while the pair (ξ, r c ) is (6, 0.76), (5.8, 0.76), (7, 0.63) for the stokeslet, stresslet, and rotlet respectively. When L = 2, M = 40, 36 and 38 for the three kernels, and the ratio M/L is kept constant as the system is scaled. We have excluded the precomputation cost in our runtimes since it is performed only once, and is easily amortized over multiple runs due to iterations or time steps when the size of the domain does not change. The figure allows us to determine the break-even value N , that is, the smallest value of N for which the Ewald summation is faster than the direct sum. In order to compare CPU runtimes, it is necessary that the simulations should use the same system and number of cores. Hence the CPU runtime study for free-space was repeated with the same parameters as above except that ξ = 7, and r c = 0.63, 0.63, 0.58 for the stokeslet, stresslet and rotlet respectively. The results are shown in Figure 5b.
The break-even values for all kernels (when excluding the precomputation cost) obtained from Figure 5 are presented in Table 4. Since the cost of the direct sum for the half-space almost quadruples in comparison to that of free-space, it benefits more from the Ewald decomposition and Spectral Ewald method and the break-even is attained much earlier. As expected, the stresslet, with the steep increase in the number of FFTs, has the largest break-even, but it is only slightly greater than the break-even value for the other kernels. These numbers underline the advantage of the Spectral Ewald method for the evaluation of the formulae (5) for the half-space.   We next study the computational run-time of different parts of the algorithm. In the left plot of Figure 6 the runtimes for real-space and Fourier-space evaluation are presented for the stresslet 3 , along with the cost of precomputation of the Green's functions. For the stresslet kernels, the precomputation involves the evaluation of the Green's function for the stresslet as well as that of the correction functions and including it in the total runtime will increase the cost by 25-33%. The usage of an optimized 3 We choose the stresslet because it is the most complicated value of the Ewald parameter ξ balances the cost of the real-space and Fourier-space evaluation and they are thus of the same order of magnitude. The right plot of Figure 6 shows a further breakdown of the Fourier-space cost into its main constituent steps, namely, Gridding, Scaling, and FFT. The scaling step is clearly the cheapest among them, and the overall results are very similar to the case of free-space despite the fact that here we perform 21 FFTs compared to only 9 for free-space.

Conclusion and further work
We have presented a fast summation method for the half-space Green's functions of Stokes flow derived by Gimbutas et al [1]. The fast summation method and its implementation follows that for free-space Green's functions by Klinteberg et al [2]. The method is based on the Ewald decomposition that recasts the sum into a sum of two exponentially decaying series: one in real-space (short-range interactions) and one in Fourier-space (long-range interactions) with the convergence of each series controlled by a common parameter. While the evaluation of the real-space component proceeded along expected lines, the presence of extra terms complicated the task for the Fourier-space component. We followed the framework of the Spectral Ewald method for free-space Stokes flow introduced recently, and exploited the structure of the terms to optimize the number of FFTs and IFFTs that need to be performed. Furthermore, we demonstrated that with very elementary modifications the truncation error estimates for free-space Stokes flow remain valid.
The implementation for the half-space does incur extra costs in comparison to the free-space in multiple ways such as the gridding of a larger computational domain, substantial increase in the number of FFTs and IFFTs that need to be evaluated but the computational savings are also greater.
Future work can take shape in one of two ways. For one, it would be beneficial to use this work in the framework of a boundary integral method for Stokes flow in a half-space motivated by a physical problem of sedimentation. A more involved, and mathematically interesting question would be to consider a 2-periodic extension with periodicity in the in-plane directions, akin to the 2-periodic extension for Stokes flow considered by Lindbo et al [13]. edges the financial support of Linné Flow Centre and thanks Ludvig af Klinteberg for helpful ideas and advice with the numerical implementation.
The self-interaction term is non-zero for the stokeslet only, given by