{'fraction_non_alphanumeric': 0.05196118075212293, 'fraction_numerical': 0.02205639083924567, 'mean_word_length': 3.9976117209387776, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast Ewald summation for Stokes flow with arbitrary periodicity

A fast and spectrally accurate Ewald summation method for the evaluation of stokeslet, stresslet and rotlet potentials of three-dimensional Stokes flow is presented. This work extends the previously developed Spectral Ewald method for Stokes flow to periodic boundary conditions in any number (three, two, one, or none) of the spatial directions, in a unified framework. The periodic potential is split into a short-range and a long-range part, where the latter is treated in Fourier space using the fast Fourier transform. A crucial component of the method is the modified kernels used to treat singular integration. We derive new modified kernels, and new improved truncation error estimates for the stokeslet and stresslet. An automated procedure for selecting parameters based on a given error tolerance is designed and tested. Analytical formulas for validation in the doubly and singly periodic cases are presented. We show that the computational time of the method scales like O(N log N) for N sources and targets, and investigate how the time depends on the error tolerance and window function, i.e. the function used to smoothly spread irregular point data to a uniform grid. The method is fastest in the fully periodic case, while the run time in the free-space case is around three times as large. Furthermore, the highest efficiency is reached when applying the method to a uniform source distribution in a primary cell with low aspect ratio. The work presented in this paper enables efficient and accurate simulations of three-dimensional Stokes flow with arbitrary periodicity using e.g. boundary integral and potential methods.


Introduction
Stokes flow, also known as creeping flow or viscous flow, is a model of fluid flow in which inertial forces are assumed to be negligible in comparison to viscous forces (i.e., the Reynolds number is very small). This is often a valid assumption for phenomena involving suspension flows on the micro-and nanoscales (often in combination with Brownian motion), such as swimming microorganisms [1], cell dynamics [2], microfluidic devices [3], gels [4,5], dynamics of nanoparticles and nanofibrils [6,7,8,9], electrolytes [10], and antibodies [11]. For Stokes flow, the Navier-Stokes equations reduce to the Stokes equations, which for an incompressible Newtonian fluid are given by −∇p(x) + µ∇ 2 u(x) + g(x) = 0, (1) ∇ · u(x) = 0. ( Here, p is the pressure, u is the fluid velocity, g is the body force per unit volume acting on the fluid, and µ is the viscosity of the fluid. We will here consider the nondimensionalized Stokes equations, which is equivalent to setting µ = 1. In boundary integral and potential methods for Stokes flow, the fundamental solution (Green's function) of the Stokes equations appears, namely the stokeslet kernel. In the three-dimensional case, the stokeslet is a 3 × 3 tensor given by where δ jl is the Kronecker delta. Also commonly used are the stresslet and rotlet kernels, which can be seen as derivatives of the stokeslet and will be introduced in section 2. In this paper, we are interested in problems with periodic boundary conditions, in which D of the three spatial directions will be periodic, and the remaining 3 − D directions will be free, in this context meaning that the domain extends to infinity and that no boundary conditions are enforced in these directions in the summation procedure. We will call the problem triply, doubly and singly periodic if D = 3, 2, 1, respectively, and free-space if D = 0. We will consider a system of N point sources of strengths f (x n ) located at positions x n . The velocity field generated by this system is given as a periodic sum over the point sources, i.e.
where the set P DP is a D-dimensional lattice containing all periodic images, to be properly defined in section 2. We assume that we want to evaluate (4) also in N target points, which may or may not be the same as the source points. The periodic sum can be computed using Ewald summation, which splits the sum into a short-range part, to be summed directly, and a smooth long-range part, to be treated in Fourier space. For the stokeslet, the split was derived by Hasimoto [12]. A decomposition parameter ξ, to be introduced in section 2, controls the Ewald summation split. For fixed ξ, computing the Ewald sums directly leads to a method with time complexity O(N 2 ). Adjusting ξ properly reduces the complexity of the direct summation to O(N 3/2 ), see e.g. [13,14]. Fast Ewald summation methods such as the Particle-Mesh-Ewald (PME) method [13] and smooth Particle-Mesh-Ewald (SPME) method [15,16] reduce the complexity further to O(N log N ), a significant improvement. The Spectral Ewald (SE) method is a PME-type method with spectral accuracy. As in all PME methods, the interactions between source points and target points are computed via a uniform grid and the fast Fourier transform (FFT). A so-called window function is used to interpolate between source/target points and the uniform grid, much like in the nonuniform fast Fourier transform (NUFFT) [17,18]. The spectral accuracy of the SE method comes from the choice of window function, which was originally a truncated Gaussian. In contrast, other methods such as e.g. SPME [15,16] use Cardinal B-splines for interpolation, which leads to algebraic accuracy.
The SE method for Stokes flow has been developed in a series of papers, first for the triply periodic stokeslet [19], doubly periodic stokeslet [20], triply periodic stresslet [21], and triply periodic rotlet [22]. The method was extended to the free-space case for all three kernels by [23]. The current paper serves to complete this development, much in the same way as was recently done for electrostatics by [24], by adding the missing pieces (singly periodic case for all three kernels, and doubly periodic case for the stresslet and rotlet), and unifying all periodic cases within a single framework. This opens up the possibility to perform efficient simulations of three-dimensional Stokes flow with arbitrary periodicity (D = 3, 2, 1, 0), using boundary integral and potential methods.
In this paper, the aperiodic directions are assumed to be free and extend to infinity. If they are instead bounded, one possibility is to explicitly discretize the boundary, and enforce boundary conditions on it, as done e.g. in [32], without modifying the underlying SE method. Another option is the general geometry Ewald-like method (GGEM) [33,34], which uses the same split into a short-range and long-range part, but treats the long-range part in real space using a mesh-based solution method. GGEM can in principle handle nontrivial boundary conditions, but involves an expensive correction solve if high accuracy is desired. A more recent solution is presented by [35], in which a Chebyshev method is used for the boundary value problem in the aperiodic direction; this method has been demonstrated for electrostatics but is expected to generalize to Stokes flow. Related work on Ewald-type methods for reduced periodicity in electrostatics includes [36,37,38]. For half-space Stokes problems, the methods by [39] or [40] can be used. The SE method has also been implemented for two-dimensional Stokes flow [41].
An alternative to Ewald-like summation methods is the fast multipole method (FMM) [42,43,44,45], which in general achieves O(N ) complexity. Unlike Ewald-type methods, which reach their highest efficiency for fully periodic problems, the FMM is most efficient and most natural to formulate in the fully aperiodic (free-space) setting. Nevertheless, the FMM has also been generalized to arbitrary periodicity [46,47,48]. An advantage of the FMM is that it is spatially adaptive, while the SE method requires a uniform spatial grid due to the FFT. Thus, the FMM will typically be faster for highly nonuniform point distributions (especially in free space), while the SE method may be faster for uniform distributions, as shown e.g. by [23,24]. Yet another alternative method is found in [49], in which the long-range interaction is represented by auxiliary sources. The principle is the same in both two and three dimensions, but it has not been demonstrated that this method would be competitive with Ewald-type methods in three dimensions.
The contribution of the current paper is, as mentioned above, to complete and unify the previous work on the SE method for the kernels of Stokes flow (stokeslet, stresslet, rotlet), by adding the singly and doubly periodic cases, and treat all periodic cases within the same framework. An important part of the unified SE method is the modified kernels that are used to treat the singular integration in cases with reduced periodicity, based on an idea by [50]. In this paper, we improve the convergence of the modified kernels used by [23] in the free-space case, and derive new modified kernels for the singly and doubly periodic cases. We also derive a new improved truncation error estimate for the stokeslet and stresslet, valid in all periodic cases, based on techniques from [23]. Analytical formulas useful for validation are derived in the singly and doubly periodic cases, completing the formulas previously derived by [20] for the doubly periodic stokeslet. The SE method presented in this paper furthermore uses the polynomial Kaiser-Bessel (PKB) window introduced by [51] and [24], and the adaptive Fourier transform (AFT) introduced by [52] for the singly and doubly periodic cases.
The paper is organized as follows. In section 2, the Ewald splits and Ewald sums are given for the three kernels. In section 3, the fast method to compute the Fourier-space Ewald sums is presented, i.e. the SE method; this section also describes the modified kernels, PKB window, and AFT. In section 4, we give error estimates for the SE method and describe an automated procedure to select the parameters of the method given an error tolerance; this procedure is also tested by numerical examples. In section 5, more numerical results follow, regarding the pointwise error, computational time and complexity, and the window function. Finally, conclusions are drawn in section 6. In the appendices, we derive analytical formulas for validation in Appendix A and Appendix B, and the improved truncation error estimate for the stokeslet and stresslet in Appendix C.

Ewald summation for Stokes flow
We consider three fundamental solutions of Stokes flow, namely the stokeslet S, rotlet Ω, and stresslet T , which are tensorial kernels given by S jl (r) = δ jl |r| + r j r l |r| 3 , Ω jl (r) = jlm r m |r| 3 , T jlm (r) = −6 r j r l r m |r| 5 , where j, l, m ∈ {1, 2, 3}. Here, δ jl denotes the Kronecker delta and jlm denotes the Levi-Civita symbol; Einstein's summation convention is used, with repeated indices implicitly summed over {1, 2, 3}. Given a point force 8πf acting on the fluid at x 0 , the Stokesian velocity field is given by u j (x) = S jl (x − x 0 )f l , the vorticity is given by ω j (x) = 2Ω jl (x − x 0 )f l , and the stress field is given by σ jl (x) = T jlm (x − x 0 )f m [53, ch. 2.2]. In boundary integral equations, the kernels S, Ω and T are multiplied by different source quantities which all give rise to velocity fields, namely u where notably the stresslet source f T lm has two indices. Commonly, the stresslet source is of the form f T lm = q l ν m , where q and ν are vectors, and this will be assumed in this paper.
We have here introduced a linear differential operator K for each kernel. To be able to treat all kernels together where appropriate, we now introduce some special notation following [23]; we will write where the kernel G may be the stokeslet S, rotlet Ω or stresslet T , and the notation G · f will be understood to mean one of S jl f l , depending on the actual kernel. (Note that G · f is always a vector quantity. ) We are interested in efficiently evaluating periodic potentials of the form generated by N sources of strengths f (x n ) located at points x n inside a box B = [0, L 1 ) × [0, L 2 ) × [0, L 3 ), called the primary cell. The set P DP of periodic images depends on the number of periodic directions D and is given by Note that the first D coordinate directions are the periodic ones; the remaining 3 − D directions are called the free directions. (The notation DP introduced here is to be read as "D-periodic".) The sum (14) may be evaluated either at one of the source locations x m (m = 1, . . . , N ), in which case the term corresponding to m = n, p = 0 is omitted, or at other arbitrary locations. In the former case, we write where the star above the second sum denotes precisely that the term p = 0 is omitted when n = m.
Since the kernels (5) decay slowly as |r| → ∞ (the stokeslet decays like 1/|r|, and the rotlet and stresslet like 1/|r| 2 ), the periodic sums (14) and (16) are absolutely convergent only in the case D = 0 (for all kernels), and D = 1 for the rotlet and stresslet (but not the stokeslet). In all other cases, the sums are either conditionally convergent (i.e. their values depend on the order of summation) or even divergent. Special care must then be taken when interpreting (14) and (16); an overview of these considerations is given by [56]. Ewald summation corresponds to a spherical order of summation.
Ewald summation is based on the idea to split the periodic sum into two parts, where the first part converges fast in real space (the "real-space part"), and the other part converges fast in Fourier space (the "Fourier-space part"). Ewald decompositions for the kernels considered here are given in section 2.1. To aid the reader, we then present the Ewald sums first in the triply periodic setting (D = 3) in section 2.2, and then for arbitrary periodicity (D = 3, 2, 1, 0) in section 2.3.

Ewald decompositions
In general, there are two different but equivalent ways to derive an Ewald decomposition, called screening and splitting [23]. The real-space part is easier to derive using the splitting approach, while the Fourierspace part is easier using the screening approach. In this paper, we use the screening approach since the Fourier-space part is our main focus.
In the screening approach, the kernel is convolved with a screening function γ(r; ξ), which is a smooth function with the property R 3 γ(r; ξ) dr = 1. The positive parameter ξ is called the Ewald decomposition parameter, and controls the decay of the resulting decomposition where G R := G * (δ − γ) is the real-space part and G F := G * γ is the Fourier-space part; here, δ is the Dirac delta distribution, and * denotes convolution. In our case where G = KA, with K being a linear differential operator and A either the harmonic H or biharmonic B, cf. (9)-(11), the decomposition can be written as Thus, we start by writing down the Ewald decompositions for the harmonic and biharmonic kernels, and may then apply the appropriate differential operator K to get the decompositions for the stokeslet, rotlet, and stresslet. For the harmonic, we select the classical Ewald screening function where the hat denotes the Fourier transform This leads to the decomposition originally derived by Ewald [57], namely where erf(·) is the error function and erfc(·) = 1 − erf(·) is the complementary error function. Here, the singular behaviour of 1/|r| is contained in H R , while the long-range behaviour is contained in H F . Note that H R decays fast as |r| → ∞ in real space, while H F is smooth and therefore its Fourier transform decays fast. By the convolution theorem, the Fourier transform (in the distributional sense) of H F is since H(k) = 4π/|k| 2 . For the biharmonic, the classical Ewald screening function does not yield a rapidly decaying real-space part, so a different screening function is needed. We will use the Hasimoto screening function [33,12] (Another option is the Beenakker screening function, see e.g. [23, Table 1], but the Hasimoto screening function yields a somewhat faster decay in both real and Fourier space, and is therefore preferred.) The Hasimoto screening leads to the decomposition By the convolution theorem, we have since B(k) = −8π/|k| 4 . Again, B R decays fast in real space, while the Fourier transform of B F decays fast. By applying the relations (9)-(11), one obtains Ewald decompositions also for the stokeslet, rotlet, and stresslet. The real-space part kernels G R become For the Fourier-space parts, we rewrite the relation G = KA, cf. (9)-(11), as G = K A, where the Fourierspace operator K is found by letting the differential operator K act on e ik·r . We have The Fourier transform of the Fourier-space part kernel then becomes G F (k; ξ) = K(k) A(k) γ(k; ξ). For the three kernels of Stokes flow, The Ewald sums are obtained by inserting G = G R + G F into (14) or (16), treating the sum associated with G R in real space and the one associated with G F in Fourier space. We do this below first in the triply periodic case, and then for arbitrary periodicity.

Triply periodic Ewald summation
Consider the case with D = 3 periodic directions. Let us first consider (14), where the evaluation point x does not coincide with any of the source locations x n . Inserting G = G R + G F into (14) yields a decomposition u 3P = u 3P,R + u 3P,F , where the real-space part is simply Note that since G R decays fast, this periodic sum is convergent. For the Fourier-space part, we use the Poisson summation formula to rewrite the periodic sum as where the box volume is |B| = L 1 L 2 L 3 and the set of discrete wavenumbers is given by In (39), the term corresponding to k = 0 (the zero mode) has been omitted since G F is singular at k = 0.
The zero mode of a Fourier series corresponds to a constant where C (0) is a constant tensor. The zero mode (41) may be added to (39), and the freedom to select C (0) reflects that the fundamental solutions are defined only up to a constant. To determine the constant, one can e.g. impose a requirement of zero mean flow [21,22], i.e.
where D j is the face of the box B = [0, For the stokeslet and rotlet, one can show [22] that (42) is satisfied when the zero mode is zero, i.e. u S,3P,F,k=0 = u Ω,3P,F,k=0 = 0. For the stresslet, however, the situation is more complicated. Setting the zero mode to zero will not result in a zero mean flow, as shown by [21]. Yet another perspective is provided by [58,Paper IV], where the Ewald decomposition for the stresslet is derived using the relation where Π j (r) = 2r j /|r| 3 is the fundamental solution for the Stokesian pressure. (The corresponding derivation in two dimensions was carried out by [59].) The resulting decomposition is the same as the one given here, i.e. (31) and (37), for k = 0. However, it is noted that by setting u S,3P,F,k=0 = 0 for the stokeslet, a mean pressure gradient will appear (in the Ewald sum for the pressure) to balance the point forces, represented by ∇ l Π (0) j (r) = 8πδ jl /|B|. This implies that the pressure itself gains a nonperiodic part which via (43) carries over to the stresslet as Thus, the stresslet decomposition will gain the additional nonperiodic term (here we abuse the notation for the zero mode to represent this extra term, which is of course not constant) Using the assumption that f T lm = q l ν m , this can be written as af Klinteberg and Tornberg [21] derived a similar term under the so-called rigid body assumption and our term (47) reduces to theirs under the rigid body assumption (48).
To summarize, the full Ewald decomposition, for an evaluation point x which does not coincide with any of the source locations, is given by where u 3P,R is given by (38), u 3P,F,k =0 is given by (39), and u 3P,F,k=0 is zero for the stokeslet and rotlet but given by (47) for the stresslet. When the evaluation point is one of the source locations x m (m = 1, . . . , N ), the term corresponding to n = m for p = 0 (the "self interaction") should be omitted, cf. (16). For the real-space part, the self interaction can be omitted directly, and the real-space part is thus given by (38) but with a star above the second summation sign (denoted by u 3P,R ). For the Fourier-space part, the self interaction must be removed explicitly, which is done by adding the term For the stokeslet, we get while for the rotlet and stresslet, u Ω,self (x m ; ξ) = u T,self (x m ; ξ) = 0. The full Ewald decomposition is then given by Here, u 3P,R is given by (38) but with a star above the second sum (to skip the term p = 0 for n = m), u 3P,F,k =0 is given by (39), u 3P,F,k=0 is zero for the stokeslet and rotlet but given by (47) for the stresslet, and u self is zero for the rotlet and stresslet but given by (51) for the stokeslet. The expressions for G R and G F for the different kernels are found in (29)- (31) and (35)- (37), respectively. In summary, (52) is used if the evaluation point coincides with one of the source locations, while (49) is used otherwise.

Ewald summation in arbitrary periodicity
Let us now consider D = 3, 2, 1, 0 periodic directions. When reducing the number of periodic directions, the real-space part is straightforward to write down, and is given by with a star above the second summation sign if x coincides with one of the source locations x m . Recall that the set P DP , given by (15), ensures that the potential is periodically summed only in the first D directions.
Since G R decays fast, the sum (53) is convergent, and it can be computed efficiently for example using a cell list, as described in section 4.1.2. The term u self , which removes self interaction from the Fourier-space part and is defined by (50), is independent of D.
In the Fourier-space part, cf. (39), the summation over discrete wavenumbers will in the free directions be replaced by integration over continuous wavenumbers (as the Fourier transform must be used instead of a discrete Fourier series). Writing down the formulas explicitly for all D (at the moment ignoring the singularity at k = 0), we have where the set of discrete wavenumbers is given by, cf. (40), In (54)-(57), we have written κ i instead of k i in the free directions to emphasize that these wavenumbers are continuous. Occasionally, we will use k P to denote the vector of discrete wavenumbers in the periodic directions, i.e. k P ∈ K D . Thus, k P means (k 1 , k 2 , k 3 ) for D = 3, (k 1 , k 2 ) for D = 2, and k 1 for D = 1. Similarly, we may use κ ∈ R 3−D to mean κ 3 for D = 2, (κ 2 , κ 3 ) for D = 1 and (κ 1 , κ 2 , κ 3 ) for D = 0. The notation k without superscript always refers to the full vector of three wavenumbers regardless of D.
We can unify (54)-(57) by introducing a "mixed" Fourier transform F DP , as follows. Let f : R 3 → C be a function that is periodic in the first D coordinate directions, and nonperiodic in the remaining 3 − D directions. We define the mixed Fourier transform of f by where the integration domain is 9 The inverse transform is given by and it holds that f (r) = F −1 DP {F DP {f }}(r). Note that this is simply a Fourier series in each periodic direction, and a Fourier transform in each free direction. We can now write (54)-(57) compactly as Let us consider the integrals that appear in the operator F −1 DP in (62) for D = 2, 1, 0. Note that G F , cf. (35)- (37), and thus the integrand, has a singularity at k = 0. This means that the integrand is singular for D = 2, 1 when k P = 0 (but not when k P = 0), as well as for D = 0 (always). In fact, the integrals may not even exist in the Lebesgue sense, but they can be interpreted as inverse Fourier transforms in the distributional sense, as seen in Appendix A.2 and Appendix B.2. For D = 2 (see Appendix A), the integrals can be evaluated analytically in both the nonsingular and singular cases (in the distributional sense in the latter case), for all three kernels. The results are of the form (with x n = (x n , y n , z n )) u 2P,F,k P =0 (x, y, z; ξ) = 1 u 2P,F,k P =0 (x, y, z; ξ) = 1 where Q 2P and Q 2P,(0) are tensors that depend on the kernel; Q 2P is given by Also for D = 1 (see Appendix B), the integrals can be evaluated analytically. The results are of the form (again with x n = (x n , y n , z n )) where Q 1P is given by (B The formulas given in (63)-(66), Appendix A, and Appendix B for D = 2, 1 are not part of the fast method to be introduced in section 3, but they serve to establish that the Fourier integrals appearing in (62) are well-defined (in the distributional sense), even in the singular case. Furthermore, they can be used to validate the fast method, by truncating the periodic sums (63) and (65) at some maximum absolute wavenumber k i,max in each periodic direction (so that only k i ∈ {2πk i /L i : −k i,max ≤k i ≤k i,max −1,k i ∈ Z} are included in the sum), and then directly sum them.
For D = 0, the integral in (57) can be seen to exist for all three kernels (also in the Lebesgue sense), for example by going to spherical coordinates. We do however not compute it analytically here, since it is typically not needed; validation is in the D = 0 case most easily done by directly summing the original sum (14).
To summarize, the full Ewald decomposition in arbitrary periodicity is given by where (68) is used if the evaluation point coincides with one of the source locations, and (67) is used otherwise. Here, the real-space part u DP,R( ) is given by (53), and the star signifies that the term p = 0 is skipped when n = m. The Fourier-space part u DP,F is given by (62) and can for D = 3, 2, 1 be further decomposed as as in (63)-(66); for D = 3, cf. (49) and note that k P = k for D = 3. Finally, the term u self in (68) is independent of periodicity; it is zero for the rotlet and stresslet, and given by (51) for the stokeslet. Again, the expressions for G R and G F for the different kernels are found in (29)-(31) and (35)-(37), respectively.

A note on zero modes, far-field behaviour and the stresslet integral identity
Let us note that the term u DP,F,k P =0 in (69) is the zero mode of the Fourier series in the periodic directions; it is a function of the coordinates in the free directions only. In the triply periodic case (D = 3), cf. section 2.2, the requirement of zero mean flow (42), was used to fix the arbitrary constant in the zero mode (41). However, zero mean flow cannot be imposed in the D = 2, 1 cases, since the computational domain is unbounded in the free directions, and, as shown (for the stokeslet) by [56], in order for the flow field (67) to be smooth it must diverge as infinity is approached in the free directions. This far-field behaviour in the free directions is controlled by the zero mode.
With the zero modes derived in Appendix A and Appendix B, the far-field behaviour of u DP,F,k P =0 in the free directions is for the stokeslet O(|z|) in D = 2, and O(log(|r|)) in D = 1. For both the rotlet and stresslet, the behaviour is O(sgn(z)) in D = 2, where sgn(·) is the sign function, and O(1/|r|) in D = 1. In D = 0, the flow field is uniquely determined by (14) and goes to zero at infinity.
In boundary integral methods, one would like the stresslet integral identity [53, eq. 2.1.12, p. 21] to hold for an arbitrary constant vector q (0) . Here, Γ is a sufficiently smooth surface enclosing the domain D encl , and D ext is the domain outside Γ; furthermore, ν is the outward-pointing unit normal of Γ, and the set P DP is as in (15). (Note that (70) upon discretization of the integral becomes precisely the stresslet flow field (14) for a constant q(y) = q (0) .) As noted by [21], the term u T,3P,F,k=0 as given by (47) is needed for the stresslet integral identity (70) to hold in the D = 3 case. We have verified numerically that (70) holds for all values of D = 3, 2, 1, 0 when the flow field is given by (67)-(69).

The Spectral Ewald method
The goal is now to compute the Fourier-space part (62) of the periodic potential in an efficient way. While it is possible to truncate and directly sum (54), (63) and (65), corresponding to D = 3, 2, 1, respectively, doing so would yield a slow method that scales at best like O(N 3/2 ), where N is both the number of sources and targets, and also with a large constant due to the evaluation of special functions in the D = 2, 1 cases. (In the D = 0 case, (14) can be summed directly, which would however scale like O(N 2 ).) To get a fast method, we instead introduce a uniform grid and compute the interaction between sources and targets via the grid, using the fast Fourier transform (FFT). This is the idea behind the class of Particle-Mesh-Ewald (PME) methods, to which the Spectral Ewald (SE) method belongs. To treat all D = 3, 2, 1, 0 within the same fast framework, we will discretize the integrals that appear in (62) for D = 2, 1, 0, rather than computing them analytically as was done in section 2.3.
The SE method, like other PME methods, has the following steps: (i) the sources f (x n ) are spread onto a uniform grid using an interpolating window function w(r), (ii) an FFT is applied on the grid, (iii) the result is scaled by the kernel G F , (iv) an inverse FFT (IFFT) is applied, and finally (v) the result is interpolated from the uniform grid to the desired target points using the window function w(r). The SE method differs from other PME methods in that the support of the window function can be varied independently of the size of the uniform grid, which allows approximation errors from the window function to be controlled separately from truncation errors from the grid, as described in section 4. The window function will be described in section 3.3; for now it can be thought of as a generic function R 3 → R with compact support.
Let us derive formulas for the steps of the SE method by introducing the window function into (62). Given a window function w(r) with Fourier transform w(k), the identity w( w) −2 w = 1 can be inserted into (62), which can then be arranged as with the mixed inverse Fourier transform F −1 DP defined as in (61). Recall that the kernel G may be the stokeslet, rotlet or stresslet, cf. (12)- (13), and the expressions for G F are given in (35)- (37). Let us define the gridding interpolant and note that the mixed Fourier transform of Φ is, by the Poisson summation formula, which appears in (71). Evaluating (72) on the uniform grid, which is efficient since w has compact support, corresponds to step (i) of the SE method, as outlined above.
Step (ii) corresponds to computing the Fourier transform F DP {Φ}(k). Let us furthermore define which corresponds to step (iii) of the SE method. Finally, by (71) and the convolution theorem, with B DP as in (60).
Step (iv) of the method corresponds to computing F −1 DP { Φ}(y) on the uniform grid. Discretizing the integral in (75) on the uniform grid leads to step (v) of the SE method.
Several integrals appear in the formulation above, namely in the operator F DP (59) of step (ii), in the operator F −1 DP (61) of step (iv) for D = 2, 1, 0, and in (75) for step (v). These will all be discretized using the trapezoidal rule and, whenever the integration domain is unbounded, truncated. The discrete sums in (61) in the periodic directions are also truncated. The operators F DP and F −1 DP can then be approximated by the FFT and IFFT, respectively.
Special care must be taken when discretizing the integrals in (61) for D = 2, 1, since the kernel G F appearing in (74) is singular for the k P = 0 mode, as mentioned in section 2.3. Also for D = 0 the integrand of (61) is singular. These singular cases are treated following [50] by introducing modified kernels, which are nonsingular in Fourier space and defined as the Fourier transform of kernels that have been truncated in real space such that they correspond exactly to the original kernels within the primary cell B containing the sources. Furthermore, for D = 2, 1 and modes k P which are close to zero, but not exactly zero, the discretization of the integrals in (61) requires upsampling, since the kernel G F varies rapidly close to κ = 0. This is handled by using an adaptive Fourier transform (AFT) introduced by [52], which uses a local upsampling factor.
Below, we first introduce the modified kernels in section 3.1; the discretized method is then presented in section 3.2. The window function is presented in section 3.3, and the AFT is described in section 3.4. In the D = 0 case, a precomputation scheme is used to accelerate computations, and this is outlined in section 3.5. Finally, the SE method is summarized in section 3.6. A large part of the method (discretization, window function, AFT) is virtually independent of the specific kernel, and is therefore the same as in [24], which treats the harmonic kernel; we here give an overview of all parts of the method, but refer to the aforementioned paper for a more detailed discussion.

Modified kernels of Stokes flow
Let us now return to the formulation (54)-(57), before the window function was introduced. We here consider the cases where the integrands are singular, i.e. the k P = 0 mode for D = 2, 1, and the integral (57) for D = 0. These integrals are We reiterate that while the integrands are singular, these integrals are well-defined in the distributional sense, as we showed in section 2.3. To treat them numerically, however, we use the idea by [50] to modify the kernel G F in order to remove the singularity at κ = 0. Recall that G F = G γ, where the screening function γ is given by (25) for the stokeslet and stresslet, and by (20) for the rotlet. Both γ and γ decay fast, and can essentially be considered to have compact support; they are furthermore nonsingular. The reason that G is singular at κ = 0 is that G decays slowly as |r| → ∞, and the idea by Vico et al. is to truncate G outside some radius R, so that its Fourier transform becomes nonsingular. The radius R is selected large enough for the truncated kernel to agree with the original kernel within the primary cell containing all sources and targets.
Finding the Fourier transform of the truncated kernel is easier if the kernel is radial. For this reason, we make use of the relation G = KA G = K A, where the radial and scalar kernel A is the biharmonic (A = B) for the stokeslet and stresslet, and the harmonic (A = H) for the rotlet; cf. (9)-(11) and (32)- (34). For each periodicity we define a truncated A according to 13 where the inverse Fourier transform is understood to be one-, two-and three-dimensional in (79), (80) and (81), respectively. Here, the rectangle function is defined by Taking the Fourier transform (with appropriate dimensionality) of A DP R , we get A DP R , and we then set Let us for clarity summarize the situation in the D = 0 case.
where we have used that Ke iκ·r = Ke iκ·r . Assuming that the support of γ is contained within a ball of radius a γ , and that |x − x n | ≤ diam(B), where diam(B) is the diameter of the primary cell B containing all sources, it can be noted that u 0P,F R agrees exactly with u 0P,F as long as R ≥ diam(B) + a γ . Similar conclusions hold for D = 2, 1. Precisely how R is selected in the SE method is described further in section 3.2.
What remains is to derive expressions for A DP R for the harmonic and biharmonic kernels for D = 0, 1, 2. Starting with the case D = 0, the kernels are three-dimensional and given by H(r) = 1/|r| for the harmonic, and B(r) = |r| for the biharmonic. The Fourier transforms of the truncated kernels were derived by [50] as where κ = |κ| = κ 2 1 + κ 2 2 + κ 2 3 . Both expressions have finite limits as κ → 0, namely H 0P R (0) = 2πR 2 and B 0P R (0) = πR 4 . As noted by [23], the truncated biharmonic (89) has a slower decay in Fourier space than the exact biharmonic B(κ) = −8π/κ 4 ; for large values of κ, we have This slower decay would affect the convergence of the SE method for kernels based on the biharmonic (i.e. the stokeslet and stresslet). (The rotlet, which is based on the harmonic, is unaffected since H 0P R (κ)/ H(κ) ∼ 1 for large κ.) One way to understand the slower decay of (89) is that the truncation of |r| at R makes the kernel discontinuous, which introduces terms proportional to 1/κ 2 in the Fourier transform. Fortunately, this can be solved in the following simple way. Let us redefine the biharmonic kernel as We have full freedom in choosing a B and b B , and may select them to get optimal decay as |κ| → ∞. From (92), we see that this corresponds to a B + R + b B R 2 = 0 and a B + 2R + 3b B R 2 = 0, which has the solution (Note that these values of a B and b B are precisely the ones that make (91) continuously differentiable everywhere outside the origin when truncated at |r| = R.) With this selection, the modified biharmonic becomes The finite limit as With (93), we have B 0P R (κ)/ B(κ) ∼ 1 for large κ. We now get the modified stokeslet and stresslet by applying (85) to A 0P R = B 0P R given by (93)-(94), and the modified rotlet by applying (85) to A 0P R = H 0P R given by (88). The arbitrary constants a B and b B that were added to the biharmonic have no effect on the stresslet flow field, but one can show, by applying K S jl = δ jl ∇ 2 − ∇ j ∇ l to the extra terms of (91), that the stokeslet S jl gains an extra term 4b B δ jl , and the stokeslet flow field (78) gains an extra contribution This is just a constant, so it can easily be adjusted for afterwards, without affecting the time complexity of the algorithm. In numerical experiments, we will subtract (95) afterwards such that the final stokeslet flow field goes to zero at infinity, as expected from the definition of the stokeslet (5).
We move on to D = 1, where the harmonic and biharmonic kernels are two-dimensional and given by [60] with ρ = r 2 2 + r 2 3 . Here, we have included some gauge constants, similar to the D = 0 case: H and B are arbitrary positive constants, while c B is an arbitrary real constant. The Fourier transforms of the truncated versions of (96)-(97) were derived by [50] for H = 1 and B = e, c B = 0. Repeating the derivations with arbitrary constants, the Fourier transforms become where κ = |κ| = κ 2 2 + κ 2 3 , and J ν (·) is the Bessel function of the first kind and order ν. We may now select H , B and c B to optimize the decay of (98)-(99) as |κ| → ∞. Noting that To get optimal decay, we should have log(R/ H ) = 0, 1 + 2 log(R/ B ) = 0, and c B − R 2 log(R/ B ) = 0, which leads to the choices With these choices, the modified kernels become with finite limits The asymptotic behaviours are H 1P which turns out to be sufficient. For the harmonic kernel, it turns out that the value of H is not important, since the error in the SE method is dominated by the k P = 0 modes in this case (nevertheless, we select H = R). For the biharmonic, the k P = 0 mode dominates the error, and selecting (102) leads to noticeable better convergence compared to e.g. the choices B = e, c B = 0 from the original derivation by Vico et al.
Before applying (84) to get the modified stokeslet, stresslet and rotlet for D = 1, let us note that the biharmonic kernel can be avoided for the first component (the one along the periodic direction) of the stokeslet flow field. The reason is that, cf. (84), and since −κ 2 B = 2 H, we get that S 11 (0, κ 2 , κ 3 ) = 2 H(0, κ 2 , κ 3 ). Thus, the first component of the stokeslet can be based on the truncated harmonic kernel H 1P R , which marginally reduces the number of floating-point operations needed to compute it. A similar result holds for the stresslet, namely which means that the first component of the stresslet flow field can also be related to the harmonic kernel. Note that the relation is a differential relation, as it must be for the equivalent of (86)-(87) to hold. We can now write down the modified kernels in the D = 1 case. The modified rotlet, which is completely based on the harmonic kernel, is given by with K Ω as in (33), H(k) = 4π/|k| 2 , and H 1P R as in (101) with the limit (103) for (κ 2 , κ 3 ) = (0, 0). The modified stokeslet is given by with K S as in (32), B(k) = −8π/|k| 4 , and B 1P R (κ 2 , κ 3 ) as in (102) with the limit (104) for (κ 2 , κ 3 ) = (0, 0). Finally, the modified stresslet is given by with K T as in (34), and where C T H and C T B are symmetric tensors with entries given by and The values of the gauge constants H , B and c B have no effect on the rotlet or stresslet flow fields. For the stokeslet, on the other hand, one can show, similar to how (B.36) was derived (see Appendix B.2), that the extra contribution to the flow field (77) is As in the D = 0 case, this is just a constant, so it can easily be adjusted for afterwards if other values of H and B are wanted in the computation of the actual flow field. In numerical experiments, we use H = R and B = R √ e in the modified kernels (to get the optimal decay), but (somewhat arbitrarily) adjust the final flow field such that (112) vanishes, by setting B = 1 and H = e −1 . (It can here be noted that the relation B / H = e is needed for the relation ∇ 2 B 1P = 2H 1P to hold.) For D = 2, the kernels are one-dimensional, and the relation (83) can in fact be simplified to the point where the biharmonic kernel can be avoided altogether. Introducing the kernels and we can write down the relations where the constant tensors D S and D Ω are given by and the constant symmetric tensor D T has entries given by Truncating the kernels (113) and (114) at |r 3 | = R and taking the Fourier transform, we get The finite limits as κ 3 → 0 are (The asymptotic behaviours are Z 2P , which turns out to be sufficient.) The modified stokeslet, rotlet and stresslet are given by respectively. Here, K S , K Ω and K T are given by (32), (33) and (34), respectively; B(k) = −8π/|k| 4 , H(k) = 4π/|k| 2 , and other variables are as above.
For the sake of completeness, we also define modified kernels in the D = 3 case, namely where A = B for the stokeslet and stresslet, and A = H for the rotlet. The stresslet has a nonzero zero mode given by (47), but that is added separately since it is in fact a nonperiodic term that depends on the target point. Since x can be moved out of the sum in (47), the zero mode can clearly be computed in O(N ) operations for N sources and targets. In summary, the modified kernels G DP R are given by (85) for D = 0 (with A 0P R given by (93) for the stokeslet and stresslet, and by (88) for the rotlet), by (107)-(109) for D = 1, by (123)-(125) for D = 2, and by (126) for D = 3. These will replace the kernel G in the expression G F = G γ in (74), i.e. step (iii) of the SE method, which will thus become

Discrete formulation
We are now ready to discretize the integrals appearing in (72)-(75). As a first step, we truncate the unbounded domain B DP (60) that appear in the integrals. The truncated domain will also be called the extended box, since it extends the box B = [0, , which contains all sources, by length δL i in each free direction. The padding of the box B by δL i in the free directions is necessary to ensure that the window function is fully contained in the extended box B DP (even when a source point is at the boundary of the box B); furthermore, δL i must be selected such that also the screening function, which is introduced through G F in (74), has decayed sufficiently at the boundary of B DP . The minimal box pad length δL i is thus determined by the window function, but depends also on the screening function; we determine its appropriate value through numerical experiments as stated in section 4.5.
The modified kernels defined in section 3.1 are required to agree with the original kernels in the extended box B DP . The truncation radius R of the modified kernels should be set as small as possible, since large values of R will make the kernels more oscillatory, cf. e.g. (88) and (93). Thus, defining the extended side length we set the truncation radius R of the modified kernels to which are the smallest possible values given that R must be direction-independent, as illustrated in Figure 1. Figure 1: Illustration of how the truncation radius R is set, two-dimensional case (corresponding to D = 1). The rectangle is the extended box B DP projected on the x 2 x 3 -plane, and R is its diagonal, given by (130). In (a), the box is a square, while in (b), it has a higher aspect ratio. Since the modified kernels are based on radial kernels, R must be the same in all directions (it must define a circle, not an ellipse). This means that the upsampling factor s 0 , given by (134), will tend to be larger for boxes with high aspect ratio. For this reason, the SE method will be most efficient when applied to a box with low aspect ratio.
A uniform Cartesian grid with grid spacing h is introduced on B DP ; the grid has M i = L i /h subintervals in each periodic direction and M i = L i /h subintervals in each free direction. It is assumed that L i and h are selected such that M i becomes an even integer in each periodic direction. In each free direction, we adjust δL i to make sure that M i also becomes an even integer, as described in section 4.5.
The uniform grid also sets the resolution and bandwidth in Fourier space. In each periodic direction, the resolution is given by and the M i discrete wavenumbers are given by wherek ∞ i := M i /2. Note that −π/h ≤ k i < π/h. In each free direction, the resolution is given by where s(k P ) is an adaptive upsampling factor depending on the periodic wavenumber vector k P , described further in section 3.4 for D = 1, 2. For D = 0, there are no periodic directions, and the upsampling factor s = s 0 in (133) is uniform and fixed to which is needed to resolve the oscillations of the modified kernel G 0P there are no free directions and thus no upsampling.) For D = 0, 1, 2, the discrete wavenumbers in each free direction are given by The upsampling factor s(k P ) is in practice adjusted upwards such that s(k P )M i becomes an even integer.
The Fourier operators F DP and F −1 DP are now approximated by the FFT and IFFT on the uniform grid, denoted by F h,DP and F −1 h,DP , respectively. For D = 1, 2, these are adaptive FFTs with adaptive upsampling factor s(k P ), described further in section 3.4. Finally, the integral in (75) is approximated by the trapezoidal rule. The steps of the SE method can now be written down in discrete from: (i) Gridding: The gridding interpolant (72) is evaluated at the grid points x j of the uniform grid: (ii) FFT: An FFT is applied, resulting in where k l is the vector of discrete wavenumbers given by (132) and (135). (iii) Scaling: Evaluate the equivalent of (127), i.e.
(iv) IFFT: Apply an IFFT, resulting in where x j are the grid points of the uniform grid. (v) Gathering: Evaluate the trapezoidal rule approximation of (75) at the target points x m , i.e.
where the sum j is over all grid points of the uniform grid.

The window function
The time has come to describe the window function w(r) appearing in the method in more detail. As already mentioned, the window function should be a function from R 3 to R with compact support; we will furthermore assume that it is given by a tensor product w(r) = w 0 (r 1 )w 0 (r 2 )w 0 (r 3 ), where the onedimensional window w 0 : R → R has compact support [−a w , a w ], with a w > 0. When evaluating the window function w(x j − x n + p) in (136), the point x j lies on the uniform grid, and the one-dimensional window w 0 should thus be evaluated in where l i is an integer such that |l i h − x n,i + p i | ≤ a w . For an example, see Figure 2, where a w = 3h. (The situation in (140) is analogous. In the special case where target points and source points are the same, it is enough to evaluate the window function for (136), since the points in (140) have the opposite sign, and we assume the window function to be even, i.e. w 0 (−r) = w 0 (r).) In general, we assume a w to be a multiple of the grid spacing h, which means that w 0 is to be evaluated in P = 2a w /h grid points; for brevity, we will call the number of evaluation points P the "window size" in the following. Clearly, the window size P will be an even integer. In the original SE method, the window function w 0 was a truncated Gaussian. Saffar Shamshirgar et al. [24] compared the truncated Gaussian window against the Kaiser-Bessel (KB) window function [62] w 0,KB (r) = where I 0 (·) is the modified Bessel function of the first kind and order 0, and β > 0 is a shape parameter.
(Note that the Fourier transform w 0,KB is used in (138).) It was found that the window size P required to achieve a given error tolerance is significantly smaller (about 40% smaller) for the KB window (142) than the truncated Gaussian window. While evaluating the KB window w 0,KB directly is expensive, it can be approximated by a piecewise polynomial with an adaptively selected polynomial degree ν, inspired by the FINUFFT library [51,63], without affecting the overall error of the SE method [24]. The polynomial approximation is constructed by interpolating the exact KB window in ν + 1 Chebyshev points in each of the P subintervals [lh, (l + 1)h] shown in Figure 2 (l = −P/2, −P/2 + 1, . . . , P/2 − 1); the piecewise polynomial is allowed to be discontinuous where two subintervals meet. The resulting approximation is referred to as the polynomial Kaiser-Bessel (PKB) window; for more details, we refer to [24]. Both the shape parameter β and the polynomial degree ν can be tied to the window size P by setting β = 2.5P, ν = min( 1 2 P + 2, 10), which means that the PKB window is uniquely determined by P . Using these parameter choices, it was shown by [24] that the PKB window is superior to the truncated Gaussian window traditionally used in the SE method, in the sense that the window size P and hence the computational time needed to achieve a given error tolerance is smaller for the PKB window. The above mentioned paper treated only the harmonic kernel, but there is good reason to believe that the conclusion holds for any kernel, and in section 5.3 we show that the PKB window is indeed superior to the truncated Gaussian also for the stokeslet, stresslet and rotlet kernels. For this reason, we will focus exclusively on the PKB window throughout sections 3 and 4.

Adaptive Fourier transform and upsampling
We now return to the adaptive upsampling factor s(k P ) introduced in (133). For D = 0, the upsampling factor is simply given by (134), while for D = 3 no upsampling is needed; what remains to explain here are the cases D = 1, 2. Upsampling is needed to increase the resolution in two cases: firstly, for the oscillatory modified kernels G DP R used for the k P = 0 mode (see section 3.1), as in the D = 0 case; secondly, for modes k P = 0 which are close to zero, where the regular kernels G have rapid variations.
Thus, in the adaptive Fourier transform (AFT) framework, the upsampling factor s(k P ) in the free directions depends on the wavenumber k P in the periodic directions, as shown in Figure 3. In practice, this means that modes k P that require upsampling in the free directions are stored separately from modes k P that do not. For simplicity, the set of periodic wavenumbers k P is partitioned into three classes, namely (i) k P = 0, (ii) the set of nonzero modes to upsample, and (iii) the set of modes not to upsample. Here, K D is given in (58), andk ∞ i = M i /2 as in section 3.2; the thresholdk * i in (144) is a parameter to be selected. The adaptive upsampling factor is where s 0 is given by (134) with the minimization over i restricted to the free directions, and s * is a parameter to be selected. Note that for a cubic box B = [0, L) 3 , (134) simplifies to s 0 ≈ 1 +  Figure 3: Schematic illustration of adaptive Fourier transform (AFT), two-dimensional example. The upsampling factor s in the free direction is different for different wavenumbers k P in the periodic direction, as given by (146). The areas of the shaded rectangles represent the number of discrete modes that are stored; for example, |K * |s * M modes are stored for the set K * .

22
The AFT is computed by first applying an FFT in the periodic directions, and then separate FFTs in the free directions for each of the three cases in (146), with the appropriate upsampling achieved by zero-padding the free directions in real space before the second round of FFTs. The inverse transform (adaptive inverse Fourier transform, AIFT) first applies three separate IFFTs in the free directions, truncates and merges the results, and finally applies an IFFT in the periodic directions. For details, we refer to [24].

Precomputation in the D = 0 case
As mentioned in section 3.2, the D = 0 case requires an upsampling factor s 0 given by (134) (s 0 ≈ 2.8 for a cubic box) in order to resolve the modified kernel G 0P R in the scaling step. This is due to the oscillatory nature of the modified harmonic (88) and biharmonic (93). However, by precomputing effective kernels corresponding to H 0P R and B 0P R , the upsampling factor in the SE method itself can be reduced to 2, which is the minimum factor needed to compute an aperiodic convolution by FFTs.
The precomputation is done as follows, where A is used to denote either the harmonic H or biharmonic B kernel: (i) The modified kernel A 0P R , i.e. (88) or (93), is evaluated on the upsampled uniform grid of size R,tr . For the biharmonic only, a "mollification" step is applied, described below. (iv) A threedimensional FFT is applied to get A 0P R,tr in Fourier space on the grid with upsampling factor 2. For more details, we refer to [23].
The result A 0P R,tr from the precomputation is stored, and the tensorial kernels (stokeslet, stresslet, rotlet) are computed by applying (85) to it. The precomputation is beneficial when computations are to be done for several different source configurations in the same computational box B, for example in a time-dependent simulation, since the precomputation is done only once at the beginning (or if B must be resized), and the lower upsampling factor of 2 can then be used in the SE method for each source configuration.
For the biharmonic, a special "mollification" step is applied to A 0P R,tr in step (iii) above. This is needed since the modified biharmonic kernel (93) decays as 1/κ 4 in Fourier space, and the truncation in step (iii) introduces sharp corners which decay as 1/κ 2 . Thus, without mollifying the result, the truncation error would increase compared to using the grid with upsampling factor s 0 throughout the SE method (i.e. without precomputation). For the harmonic (88), and indeed for the more slowly decaying biharmonic (89) used by [23], there is no difference, since those kernels already decay as 1/κ 2 before truncating. The mollification for the biharmonic (93) is done in a tensor product fashion by introducing a (4 + 4)-point transition band furthest away from the origin in each spatial direction, and multiplying A 0P R,tr by a smooth mollifier in the transition band. More specifically, we compute where indices range over j = 1, . . . , 2M 1 ; l = 1, . . . , 2M 2 ; m = 1, . . . , 2M 3 , and µ (ν) are vectors given by This assumes that the origin is at index M ν + 1 in each direction. The mollifier function f µ is selected as a sum of Gaussians f µ (t) = e −a 2 t 2 + e −a 2 (t−7) 2 , with a = 2 7 log(10 2 ), based on numerical experiments; this function improves the Fourier-space decay of A 0P R,tr,mollified enough to not introduce an increase in the truncation error (cf. Figure 4 in section 4.1). While it may seem like selecting the transition band from the original 2M ν points, rather than extending the grid, would increase the risk of polluting the result, this does not appear to be a problem in practice, as shown by the success of the parameter selection procedure presented in section 4.

Summary of the SE method
The Spectral Ewald algorithm for computing the Fourier-space part (62) of the periodic potential (14) at arbitrary target locations x m ∈ B (which may or may not coincide with source locations) can now be summarized as in Algorithm 1. As usual, the arbitrary kernel G denotes one of the stokeslet S, rotlet Ω and stresslet T , cf. (13), while A denotes the harmonic H (for the rotlet) or biharmonic B (for the stokeslet or stresslet). For D = 0, the precomputation step described in section 3.5 is done separately, in advance. In particular, for the stokeslet and rotlet kernels, the source strengths f (x n ) are vectors in R 3 , and the gridding step is done componentwise such that Φ(x j ) is also in R 3 . Since the output u DP,F h (x m ; ξ) is also a vector in R 3 , a total of 3 three-dimensional FFTs and IFFTs are needed (one 3D FFT+IFFT for each vector component).
For the stresslet kernel, the source strengths f (x n ) are 3 × 3-tensors, and gridding is again done componentwise; in the scaling step, a 3 × 3-tensor is transformed into a vector in R 3 , such that the output u DP,F h (x m ; ξ) is a vector. Thus, the stresslet requires 9 three-dimensional FFTs and 3 three-dimensional IFFTs.

Error estimates and parameter selection
The Spectral Ewald (SE) method as stated in Algorithm 1 has a multitude of parameters, such as the decomposition parameter ξ, uniform grid spacing h, box paddings δL i , window size P , and several upsampling parameters. There is an inherent freedom in the choice of ξ, which serves to shift computational effort between the Fourier-space and real-space parts. The optimal value of ξ, which minimizes the total computation time, depends on the implementation and machine, and must be determined by experiment.
This section serves to establish an automated procedure to, given ξ, select all other parameters, such that the error of the method is as close as possible to a given error tolerance (to be made more precise below). The error will be measured in the root mean squared (rms) sense; we define the absolute rms error as where N t is the number of evaluation points x m , and u h is the approximation of the potential given by the SE method, while u ref is a reference potential with negligible error. (The nature of the reference potential will be stated in each case below.) We also define the relative rms error as In (150) and (151), the potential u may be either the full periodic potential u DP given by (14), the Fourierspace part u DP,F given by (62), or the real-space part u DP,R given by (53); to clarify what is meant in each case we will write E F rms for the Fourier-space part error and E R rms for the real-space part error (and similarly E F rms,rel and E R rms,rel for the relative errors). The automated parameter selection procedure to be described in this section will be given an absolute error tolerance τ rms , and is to select parameters such that the actual error is within one order of magnitude of the given tolerance, i.e. τ rms 10 ≤ E rms ≤ 10τ rms .
To simplify formulas, we will in sections 4 and 5 restrict ourselves to the case where the primary cell B is a cube of side length L, i.e. L 1 = L 2 = L 3 = L. Note that this is not a restriction of the method itself, nor of the implementation, which works for any rectangular cuboid B. Errors in the SE method come from different sources, such as truncation errors caused by truncating the Ewald sum at some maximum wavenumber k ∞ = π/h related to the grid spacing h (cf. section 3.2), and approximation errors caused by approximating integrals by discrete sums (both in (140) and in the Fourier transforms in free directions) and discretizing the window function (cf. section 3.3). Due to e.g. aliasing, these errors are interdependent, and the total error E rms is in general not equal to the sum of the individual error sources considered in isolation from each other. We will refer to this interdependence as "error pollution", and it needs to be taken into account when selecting parameters.
We will start by motivating the different aspects of the parameter selection procedure: estimates for truncation errors, the potential rms value, and approximation errors, are given in sections 4.1, 4.2, and 4.3, respectively; error pollution and its implications are discussed in section 4.4; the selection of the box padding and upsampling parameters in the free directions is described in section 4.5. Finally, the parameter selection procedure is summarized and demonstrated in section 4.6.

Truncation error estimates 4.1.1. Fourier-space part truncation error estimates
The nonzero grid spacing h used for the uniform grid in the SE method corresponds to truncation of modes with absolute wavenumbers greater than k ∞ := π/h in the Fourier-space Ewald sum (62). For the harmonic kernel, there exists an excellent estimate by [64] for the error caused by omitting the modes above k ∞ (assuming other error sources are negligible, such as when using direct summation of the Ewald sum). For the rotlet, an equally excellent estimate was derived by [22], namely where f (x n ) are the rotlet source strengths. Here, E F rms,trunc is understood to mean the absolute rms error in the Fourier-space part potential caused by truncation (i.e. when other error sources are negligible). The estimate (153) works well for any periodicity D = 0, 1, 2, 3, just like the harmonic estimate does.
For the stokeslet, the situation is a bit more complicated. Since the modified stokeslet (used for D = 0 and for the k P = 0 mode for D = 1, 2) is based on the truncated biharmonic kernel, special care is needed to ensure that the decay is fast enough as κ → ∞, as described in section 3.1. If this is not done, the truncation error becomes much larger in the D = 0 case, as illustrated by the black crosses in Figure 4; a truncation error estimate specifically for this case was derived by [23]. However, using our optimized modified kernel from section 3.1, the truncation error becomes independent of D, as it is for the rotlet. Truncation error estimates for D = 3, 2 have previously been derived by [19] and [20]. In Appendix C, we derive an improved, sharper estimate using the technique by [23] adapted to D = 3; since the truncation error is now independent of D, this estimate is valid for any periodicity. Thus, with Q given by (154) where f (x n ) are the stokeslet source strengths. The situation for the stresslet, which is also based on the biharmonic, is very similar to the stokeslet, and in fact the stresslet estimate can be related to the stokeslet estimate. Previously, an estimate for D = 3 has been constructed by [21] using curve fitting, and an estimate for D = 0 was derived by [23]. We derive an improved estimate in Appendix C, which is valid for any periodicity. Our truncation error estimate is where f (x n ) are the stresslet source strengths. The estimates are illustrated by an example with random sources in Figure 4. Whenever we refer to random sources in this paper, we will, unless otherwise stated, mean that the locations of the sources x n are uniformly distributed within the primary cell B, and that each component of the source strengths f (x n ) is uniformly distributed in the interval [−a, a], where a is adjusted a posteriori to get the value of Q stated in the example.
Note that the estimates (153)-(156) can all be solved for k ∞ = π/h, in most cases using the Lambert W function, defined as the solution to W (t)e W (t) = t; thus, k ∞ , and hence the grid spacing h, can be computed given ξ and an absolute error tolerance. Once the desired grid spacing h = π/k ∞ has been computed from the estimates, the grid size M , i.e. the number of subintervals of the uniform grid in each periodic direction, can be computed as M = L/h. However, M must at the very least be an integer, and we will require it to be a multiple of f M , where f M is some positive integer. Thus, the grid spacing is adjusted according to where h actual is the grid spacing used in the SE method, h target is the one given by the truncation error estimate, and · is the ceiling function. Typically, f M is set to a small power of two to increase the efficiency of the FFTs. Throughout section 4, we use f M = 2 in the numerical examples.

Real-space part truncation error estimates
To efficiently compute the real-space part potential (53), a cut-off radius r c > 0 is introduced, and only terms for which |x − x n + p| < r c are included in the sum. The primary cell B is divided into a uniform Cartesian mesh of rectangular cuboid subcells, such that the side lengths of each subcell are no less than r c . A list of the points in each subcell is constructed, so that the computation of (53) can be done efficiently, only considering the 27 neighbouring subcells of each evaluation point. Truncating the real-space Ewald sum at r c introduces an error, which is well described by existing estimates for the stokeslet [20], stresslet [23], and rotlet [22]: where Q is given by (154) for the stokeslet and rotlet, and by (157) for the stresslet. Estimates (160)-(162) work well for any periodicity D = 0, 1, 2, 3, for all kernels. The estimates are illustrated in Figure 5. They can all be solved for r c , again using the Lambert W function, so that they allow us to determine r c given ξ and an absolute error tolerance.

Potential rms value estimates
The truncation error estimates in section 4.1 are for the absolute error (150), and this is the reason that the parameter selection procedure takes an absolute error tolerance, rather than a relative tolerance. Of course, once the potential has been computed, one can compute the relative error (151) a posteriori (dividing the absolute error by the rms value of the potential), and relate absolute and relative errors. If one can find an estimate of the rms value of the potential, i.e. the denominator of (151), it would be possible to make an a priori estimate of the relative error, such that a relative tolerance can be directly given.
For the Fourier-space part of the potential, it turns out that one can estimate the rms value using only the box side length L, source strength quantity Q and decomposition parameter ξ, as we outline below. Not only can this be used to estimate the relative error of the Fourier-space part, but it will also turn out to play a crucial role in the error estimates to follow in sections 4.3 and 4.5, and is thus a critical piece of the parameter selection procedure. In contrast, the full potential or real-space part cannot have such a simple rms value estimate, since they depend heavily on how close to the source points the potential is evaluated, and thus have large variations in space. While it might be possible to derive more complicated estimates for the full or real-space part potentials, for example taking into account the distance to the source points, having such estimates is not critical for the parameter selection procedure, and we will therefore not do it here.
Shifting our attention to the Fourier-space part, we want to find, for each kernel, a quantity U (L, Q, ξ) such that Using numerical experiments and curve fitting, we have determined the estimates with Q as in section 4.1. Here, the functions f U for each kernel are given by (again from curve fitting) The constants C U in (164)-(166) depend on the specific configuration of source points and evaluation points.
Assuming that the points are drawn from a uniformly random distribution (as described in the paragraph following (157)), the constants C U have approximately symmetrical distributions, with mean and standard deviation given in Table 1. The data in this table is generated using 160 different random point systems per kernel, and using ξL = 5, 10, 15, 20, 25, 30 for each system; furthermore, L = 1 and Q = 1 since C U does not depend on these parameters. This was done in the fully periodic (D = 3) case, but the order of magnitude of the Fourier-space part potential does not depend strongly on the periodicity, so (164)-(166) can be used in any periodicity. When using (164)-(166) in the parameter selection procedure, the mean value of C U from Table 1 is used. These potential rms value estimates allow us to approximate the relative error of the Fourier-space part as E F rms,rel ≈ E F rms /U (L, Q, ξ). We do not account for further tests of these estimates here, but they will be implicitly tested as part of the parameter selection procedure in section 4.6.

Approximation error estimates
We will now characterize approximation errors introduced by the window function, and will therefore for the moment assume that the truncation errors described in section 4.1, and any errors introduced by additional approximations in the free directions, are negligible. In [24], the (absolute) window approximation error estimate E F rms,window ≈ 5U (L, Q, ξ) e −2πP 2 /β + erfc( β) was constructed for the exact KB window (142). Here, P is the window size and β is the shape parameter; furthermore, U (L, Q, ξ) is an estimate of the rms value of the Fourier-space part potential, cf. section 4.2 and in particular (164)-(166). In (170), the exponential term approximates the error caused by discretizing the window function and applying the trapezoidal rule to it in (140); the erfc term approximates the error caused by the discontinuity of the KB window (142) at |r| = a w = hP/2. It was shown that the two terms are approximately balanced if β ≈ √ 2πP ≈ 2.5P , at which point (170) simplifies to This window approximation error is independent of both the kernel and the periodicity, assuming that other error sources are negligible. The estimate (171) is illustrated in Figure 6 (a), and is in excellent agreement with the actual error. In our parameter selection procedure, we can use (171) to compute the window size P given ξ and an error tolerance.
In the SE method, we do not use the exact KB window, but a piecewise polynomial approximation of it, the PKB window (see section 3.3). In Figure 6 (b), we show that when the polynomial degree ν is selected according to (143), the PKB window introduces no further error compared to the exact KB window.

Error pollution
So far, we have considered errors, such as the truncation error and window approximation error, in isolation. It turns out that this is not sufficient, as shown by Figure 7. This figure has two sets of error curves: (I), for which the grid size M is selected large enough to make the truncation error negligible (so that the window approximation error dominates), and (II), for which M is selected to make the truncation error approximately of the same size as the window approximation error for each given value of P (according or M = 48 (stresslet), corresponding to a relative truncation error below 10 −23 (negligible). Note that the quantity on the y-axis is an approximation of the relative window approximation error E F rms,window,rel . The reference potential has the same parameters but is computed using the exact KB window with P = 20.
to the error estimates from sections 4.1-4.3). It is seen for (II) that the actual error is significantly larger than one would expect from the error estimates (by simply summing the individual error sources). From Figure 7 (a), it is clear that the convergence is still spectral in P for (II), but with a slower rate than expected from the window approximation error estimate (171). Thus, if both M and P are set from error estimates, the error tolerance will not be met. We call this interdependence of the truncation and approximation errors "error pollution".
We will not attempt an exhaustive analysis of the mechanics behind error pollution here. However, we note that it is related to aliasing errors that appear when evaluating the discrete Fourier transform coefficients of the window function. Aliasing causes the actual approximation error to depend not only on P , but also on M , an effect that is responsible for the "error pollution" and is missing from the approximation error estimates found in section 4.3. A more detailed analysis of the approximation error would include the effect of approximating the continuous Fourier transform of the window function by the discrete Fourier transform, which would take aliasing, and thus the "error pollution", into account; this is left as future work. A related analysis of aliasing errors in the nonuniform fast Fourier transform (NUFFT), although not directly applicable to our situation, is found in [65,66].
Error pollution can be reduced either by increasing M (which reduces aliasing errors) or P (which reduces the approximation error, albeit slower than expected from the error estimates). In practice, the best strategy with respect to minimizing the run time of the method is probably to increase both M and P moderately. We have found heuristically that if one either increases M to 1.05M and P to P + 4, or M to 1.1M and P to P + 2, compared to the values given by the error estimates, the total error decreases to the level one would expect from the estimates. Which of these rules is the fastest will depend on the periodicity D, since, as will be seen in section 4.5, increasing P also increases the extended grid size M in the free directions, thus increasing the cost of FFTs and the scaling step. It is thus more expensive to increase P the smaller D is. We have found that a suitable adjustment rule to mitigate error pollution is (h target , P target ) = (h estimate /1.1, P estimate + 2), for D = 0, (h estimate /1.05, P estimate + 4), for D = 1, 2, 3, where, as usual, h = L/M . To reiterate, the dependence on D is only to try to minimize the overall run time for each D. The rule (172) is used whenever both h and P are selected from the error estimates, unless otherwise stated. It is verified together with the rest of the parameter selection procedure in section 4.6.

Box padding and upsampling
The estimates covered in sections 4.1-4.4 are enough to select all parameters that appear in the fully periodic (D = 3) case. When some, or all, spatial directions are free, i.e. in the D = 2, 1, 0 cases, additional parameters appear in the method: the box paddings δL i , upsampling factors s 0 and s * , and upsampling thresholdk * i . Since we assume here that the primary cell is a cube, the parameters will be the same in all directions, i.e. δL i = δL andk * i =k * , i = 1, 2, 3. The parameters that we discuss here are to some extent interdependent, since, as seen by (133), both the upsampling factor s (i.e. s 0 or s * ) and the box padding δL (which appears in the side length L = L + δL of the extended box) will influence the wavenumber resolution in the free directions. Thus, if δL is selected larger, s can be selected smaller, and vice versa. Below, we first describe how δL is selected, and then how the upsampling parameters s 0 , s * andk * are selected, which is also the order that the parameters are selected in the final parameter selection algorithm.
The box padding δL relates the side lengths of the original box L and the extended box L = L + δL, in the free directions. Recall from section 3.2 that the reason that the original box must be padded is that (i) the window function must be fully contained in the extended box, and (ii) the screening function must have decayed sufficiently at the boundary of the extended box. Point (i) is guaranteed by requiring that δL ≥ 2a w = hP , where a w is the halfwidth of the window function (see Figure 2). Point (ii) is less straightforward, since the screening function does not truly have compact support, so it must be truncated at some level where it has decayed "sufficiently". We have found it useful to let δL ≈ λhP , where the factor λ ≥ 1 depends on the kernel (and thus on the screening function); this naturally connects points (i) and (ii), and furthermore links δL to the window size P , and thus to the overall precision of the method (so that δL becomes larger when more precision is required). Since the Hasimoto screening function (25) decays slower than the Ewald screening function (20), λ is expected to be larger for the stokeslet and stresslet, than for the rotlet.
In practice, we want the extended grid size M = L /h to be a multiple of f M , just like M is (cf. (158)). Therefore, we first set M , and then compute L = hM (and hence δL = L − L), rather than the other way around. The extended grid size is set to where f M is some positive integer and · is the ceiling function; here, θ is a safety threshold and θ = 8 for the stokeslet and stresslet, while θ = 0 for the rotlet. Note that (173) is approximately the same as setting δL/h ≈ λP . It remains to select λ, but we will postpone that for a few paragraphs until the upsampling parameters have been discussed. Recall from section 3.4 that upsampling is needed both to resolve the modified kernels used for the k P = 0 mode, which is upsampled by the factor s 0 , and to resolve the regular kernels for modes k P ∈ K * , given by (144), which are upsampled by the factor s * . The parameterk * controls the size of the set K * . It has been established by [23] that s 0 should be selected according to (134), where i in the minimization is restricted to the free directions (D + 1 ≤ i ≤ 3). In practice, we round s 0 upwards to one decimal, i.e. we use 10s 0 /10.
Ref. [24] determined a rule for selecting s * andk * for the harmonic kernel, namely where E F rms is the absolute rms error and U (L, Q, ξ) is an estimate of the rms value of the Fourier-space part potential. It turns out that (174)-(175) work well also for the stokeslet, stresslet and rotlet, using the appropriate U (L, Q, ξ) for each kernel according to (164)-(166). Finally, we want the upsampled grid sizes s 0 M and s * M to be multiples of f M , and will therefore adjust both s 0 and s * upwards such that this holds.
Note that the formulas (134) and (174)-(175) to select s 0 , s * andk * depend on the extended box size L and grid size M , which makes the upsampling parameters depend on the box padding δL. The only missing piece of our parameter selection algorithm is the value of λ which is needed to compute M , and thus δL, from (173). We determine λ manually by numerical experiments, by considering a plot such as Figure 6 (a), initializing λ to 1 (which makes the errors large) and then gradually increasing λ until the errors follow the estimates as in Figure 6. At every point of this process, the upsampling parameters are selected according to (134) and (174)-(175).
The resulting values of λ are given in Table 2. For D = 0, the smallest value of λ given by the process just described is used; as expected, λ is larger for the stokeslet and stresslet than for the rotlet, due to the different screening functions used for these kernels. (The stokeslet and stresslet also have slightly different λ since the relation between ξ, h and P is different for the different kernels.) For D = 1, 2, we select a larger λ than strictly necessary, namely λ = 2.4 for all kernels. This allows us to use slightly smaller values for s * andk * (since they are inversely related to M and thus to λ, cf. (174)-(175) and (173)), which reduces the overall computation time of the method slightly.  4.6. Summary of the parameter selection procedure For a given kernel (stokeslet, stresslet or rotlet) and periodicity D (0, 1, 2 or 3), the automated parameter selection procedure can be summarized as follows. The input consists of the primary cell B = [0, L) 3 , the source strength quantity Q as given by (154) or (157), the decomposition parameter ξ, and an absolute error tolerance τ rms . The general structure of the procedure is the same as in [24], with the following steps: 1. Given L, Q and ξ, compute the real-space part cut-off radius r c from estimate (160), (161) or (162), depending on the kernel, with E R rms,trunc = τ rms . (This is the only step needed for the real-space part.) 2. Compute a preliminary grid spacing h estimate = π/k ∞ by computing k ∞ from one of the estimates (153), (155) or (156), depending on the kernel, with E F rms,trunc = τ rms . 3. Compute a preliminary window size P estimate from the estimate (171), with E F rms,window = τ rms , and with U (L, Q, ξ) for the kernel in question given by (164)-(166); the constant C U is given by its mean value from Table 1. 4. Apply the error pollution adjustment rule (172) to adjust the preliminary values h estimate and P estimate , yielding h target and P target . Then apply rounding according to (158), with f M set to a small power of two, to compute the actual grid spacing h actual . This also yields the grid size M actual = L/h actual in each periodic direction. 5. Compute P actual = 2 P target /2 , then compute the shape parameter β and polynomial degree ν from (143). If D = 3, stop here; the parameter selection is complete. 6. Compute the extended grid size M in the free directions from (173), with λ given by Table 2, and using f M from step 4. Compute L = hM . This also sets δL = L − L and the truncation radius R through (130). 7. Compute the upsampling parameters s 0 from (134) with D + 1 ≤ i ≤ 3, s * from (174), andk * from (175). Here, E F rms = τ rms , and U (L, Q, ξ) is as in step 3 above. Round s 0 to 10s 0 /10, then adjust both upsampling factors such that s 0 M and s * M are multiples of f M .
The grid sizes in steps 4, 6 and 7 are consistently rounded upwards to multiples of f M , with f M set to a small power of two in order to increase the efficiency of FFTs. Typically, a small value for f M (such as 2) will put the error of the SE method close to the target error tolerance, while a slightly larger value (such as 4 or 8) may lead to more consistent run times but give a smaller error than expected. Here, in section 4, we set f M = 2 in the numerical examples, while in section 5 we will set f M = 4. (The default value in the SE package is f M = 4, but the optimal value may depend on the FFT implementation.) To verify that the parameter selection procedure works, we will consider some examples. Here, we consider only the error in the Fourier-space part of the potential; the full potential will be considered in section 5. First, we consider a single particle system per kernel, consisting of N = 1000 sources in a box of side length L = 1; the sources are also the evaluation points. The sources are random and selected such that Q = 1. We pick ξ = 10 and set the absolute tolerance τ rms to different values between 10 −1 and 10 −15 . For each kernel and periodicity, the actual Fourier-space part error is plotted versus the tolerance in Figure 8 (a)-(c). We show that the error pollution adjustment (172) is needed by plotting data both when it is used (colored unfilled symbols) and when it is not used (light gray filled symbols); the error in the latter case is too large. With the error pollution adjustment, the error is within one order of magnitude of the tolerance, as desired, until it flattens out.
Using the potential rms value estimate U (L, Q, ξ), it is also possible to prescribe a relative error tolerance τ rms,rel for the Fourier-space part, using the relation τ rms = U (L, Q, ξ)τ rms,rel . In Figure 8 (d)-(f), we plot the relative error versus the relative tolerance. Again, the error is within one order of magnitude of the tolerance when the error pollution adjustment is used, which shows that U (L, Q, ξ) is a good estimate of the potential rms value.
To further test the robustness of the procedure, we consider 90 different particle systems per kernel, given by all combinations of L = 0.1, 1, 10; Q = 0.01, 1, 100; and 10 different random seeds. For each system, we set ξ such that ξL = 10, 20, 30 and compute the actual Fourier-space part error for a few different tolerances (again evaluating the potential at the source points). The result is shown in Figure 9. In the vast majority of cases, the error is within one order of magnitude of the tolerance, the only exception being that the error in the D = 3 case sometimes becomes "too small", i.e. smaller than one tenth of the tolerance. This is of course not a problem, so we conclude that the parameter selection procedure works well. Other parameters than those already mentioned are set by the parameter selection procedure. The reference potential is computed using the same method with a relative error tolerance of 10 −17 . The rms value of the potential varies with L, Q and ξ, and ranges from 0.02 to 1000 for the stokeslet, from 0.006 to 10 000 for the rotlet, and from 0.02 to 40 000 for the stresslet.

Numerical results
Having established an automated procedure for selecting the parameters of the SE method, we move on to other numerical results. We start by considering the pointwise error in section 5.1, then measure the computational time of the method in section 5.2, and finally in section 5.3 compare the PKB window function with the truncated Gaussian window that appeared in previous iterations of the SE method.
All computations are done on 6 cores on a machine with an Intel Core i7-8700 processor, running at 4.6 GHz with 32 GB of memory. The SE method is implemented in Matlab, with critical routines written in C and called through Matlab's MEX interface. Fourier transforms are computed using Matlab's fft function, which is based on the FFTW library [67]. Throughout section 5, we let the grid sizes be multiples of f M = 4, since we have observed that this leads to more consistent timings than f M = 2. The reason is that the FFT can be substantially faster when the grid size is a multiple of 4.  [68] showed that the pointwise error has a weak dependence on the coordinate in the free direction. The same phenomenon can be seen for the kernels considered here, as shown in Figure 10 for the stokeslet (Fourierspace part). For D = 3, the error has the same size in the whole primary cell, but for D = 2, 1, the error is slightly larger at the edges of the cell in the free directions, than at the middle of the cell. However, the variation is not large, and the rms error is still within one order of magnitude of the tolerance (as already shown in section 4.6).

Timing results
We will here study how the computational time of the method as a whole (Fourier-space and real-space parts, i.e. the full potential) varies with the number of points N in the primary cell and the error tolerance τ rms . We fix the box side length L = 1 and source strength quantity Q = 1, since these do not affect the computational time. The potential is evaluated at the source points, i.e. N t = N . The precomputation step for the D = 0 case, described in section 3.5, is excluded from the timing since it can be reused for multiple solves.
First, we fix τ rms = 10 −8 and vary the number of points N , i.e. we measure the time complexity of the method. This is done for two different types of random systems, shown in Figure 11: the uniform system that we have already used above, and a clustered system consisting of denser point clouds (with points only present in the subdomains [0, 1/3) 3 and [2/3, 1) 3 ). As N increases, the point concentration N/L 3 increases, and the cut-off radius r c must be adjusted to prevent the real-space part from scaling like O(N 2 ). We set r c such that the expected average number of points within a ball of radius r c , i.e. N rc := 4 3 πr 3 c N/L 3 , is kept fixed as N increases. This way, the time of the real-space part will scale like O(N ). Given r c , we determine ξ from the truncation error estimates (160)-(162), with E R rms,trunc = τ rms = 10 −8 , and all remaining parameters are then set according to the parameter selection procedure. Roughly, we will have ξ ∼ r −1 c ∼ N 1/3 , and furthermore the grid size M ∼ k ∞ ∼ ξ. The gridding and gathering steps scale like O(N ), while the FFTs scale like O(M 3 log M 3 ), and the scaling step like O(M 3 ). Since M ∼ N 1/3 , the Fourier-space part, and therefore the method as a whole, will scale like O (N log N ).
The expected average number N rc of particles within a ball of radius r c (from which ξ is computed) is set to balance the real-space part and Fourier-space part run times, so that the total computation time is close to minimal. The optimal value depends on the particle system (uniform or clustered), as well as the periodicity and kernel (and the computer), and is given (approximately) in Table 3. Note that N rc is simply the volume of an r c -ball times the average particle concentration N/L 3 in the whole primary cell (including also the empty regions for the clustered system; N/L 3 has the same value for the uniform and clustered systems). The fact that N rc is selected smaller for the clustered system than for the uniform one reflects that r c should be smaller, since the particle clouds are denser. The reason that N rc is selected larger for lower periodicities D is that the Fourier-space part becomes more expensive, so more computational effort is shifted into the real-space part.  Rotlet  2500  950  450  400  800  300  170  120  Stresslet  6000  1600  800  800  3000  1200  370  300 The computation time when increasing the number of points N , while keeping N rc fixed, is shown in Figure 12 for the uniform system, and in Figure 13 for the clustered system. It can be noted that the computation time follows the expected scaling O(N log N ) for both systems, and that the computation time is a bit larger for the clustered system than for the uniform one. The uniform particle distribution is the optimal one for the SE method, since there is no spatial adaptivity in the Fourier-space part method.

Uniform system Clustered system
As a second test, we fix N = 10 5 and measure the computation time when varying the error tolerance τ rms between 10 −2 and 10 −14 . In this test, we use a fixed decomposition parameter ξ, corresponding to the value selected (from N rc ) in the time complexity test for N = 10 5 and τ rms = 10 −8 ; the value of ξ is given per kernel, periodicity and particle system in Table 4. Keeping ξ fixed as the tolerance τ rms is varied allows the run time of both the real-space and Fourier-space parts to be adjusted (through r c and M ), so that the computational balance between the two parts is kept approximately the same. Parameters are set according to the parameter selection procedure, as in the complexity test.
The result is shown in Figure 14 for the uniform system. (The clustered system has a very similar behaviour, but with a factor 3 to 4 larger computational times, as already seen in the complexity test. It is therefore omitted here.) For each kernel, it can be seen that the computation time for D = 2 is very similar to D = 3, while the time for D = 1 is at most twice as large; furthermore, the time for D = 0 is around three times as large as D = 3. It can also be seen that the stresslet requires more time than the other kernels (up to twice as much) for the same relative error; this is because the stresslet source strengths have more components, which for example leads to more expensive gridding and more FFTs.     Table 3 for N = 10 5 and τrms = 10 −8 .

Window function comparison
Finally, we compare the PKB window function used in this paper with the truncated Gaussian window traditionally used in the SE method. The truncated Gaussian (TG) is given by (cf. (142)) where a w = hP/2, and α is a shape parameter which has the value α = 0.91(π/2)P . In the scaling step, the Fourier transform of the window function is needed. Rather than using the exact Fourier transform w 0,TG of the truncated Gaussian, which involves the error function erf(·), we prefer to use the Fourier transform of the untruncated Gaussian w 0,UG (r) = e −α(r/aw) 2 , i.e.
The reason is that w 0,UG is cheaper to evaluate than w 0,TG , and the difference between them is guaranteed to be around the selected error level due to how parameters are selected. Thus, we use w 0,TG in the gridding and gathering steps, but w 0,UG in the scaling step.
To compare the PKB and TG windows, we consider for each kernel a system of N = 10 5 uniformly random sources, with Q = 1, in a cubic box of side length L = 1. The parameter ξ is selected as in Table 4 (uniform system). The parts of the SE method that depend most strongly on the window function are the gridding and gathering steps, which have run time O(P 3 ) with respect to the window size P . Ref. [24] showed, for the harmonic kernel, that P can be reduced about 40 % for the PKB window compared to the TG window, for the same error. Thus, the gridding and gathering time will be smaller for the PKB window than for the TG window. In Figure 15 we show that this is indeed true also for the three kernels that we consider here: The gridding and gathering time for the PKB window is significantly reduced compared to the TG window, since the window size P can be selected smaller. This effect becomes more pronounced for stricter error tolerances. The stokeslet and rotlet have identical gridding and gathering steps; the only reason that the run time is larger for the rotlet than for the stokeslet in Figure 15 is that ξ, and therefore M , is selected larger for the rotlet, in accordance with Table 4.
Here, we have focused on the run time of the gridding and gathering steps of the SE method. The FFT and IFFT steps do not depend on the window function, while the scaling step has only a weak dependence on the window (since the Fourier transform w 0 must be evaluated). Thus, the PKB window will typically lead  Table 4. Only the fully periodic (D = 3) case is included here, but gridding and gathering times do not depend strongly on periodicity. The grid size is M = 196 for the stokeslet, M = 192 for the stresslet, and M = 208 for the rotlet, which corresponds to a negligible absolute truncation error around 10 −20 . The data points are obtained by letting the window size P vary from 2 to 14 for the PKB window, and from 2 to 24 for the TG window (right to left), computing the error and the time for each P . Each data point is the average of 10 different point configurations. The rms value of the Fourier-space part of the potential is around 2 for the stokeslet, 50 for the stresslet, and 20 for the rotlet. The reference potential is computed using the same method, with the PKB window and P = 20.
to smaller computation times than the TG window, also when the total computation time of the method is considered. The percentual total time gain of using the PKB window will depend on how much time is spent in the gridding and gathering steps compared to the other steps, which depends on different factors such as the number of points N , decomposition parameter ξ, periodicity D, and the error tolerance τ rms (which sets M and P ). In the example considered here, the gridding+gathering time makes up from 40 % (for P = 2) to 70 % (for P = 24) of the total Fourier-space part run time.
We would like to point out that the SE method has been shown to be competitive with other fast summation methods, such as the fast multipole method (FMM); for instance, this has been done for the harmonic kernel by [24], and for the stokeslet by [23]. In the latter case, the comparison between SE and FMM was done in free space (D = 0), which is the most expensive periodic case for the SE method, but the cheapest for the FMM. The comparison was done with the TG window function in the SE method; as we have shown here, the PKB window can further reduce the run time of the SE method.

Conclusions
In this paper, we have presented the Spectral Ewald (SE) method, a fast Ewald summation method, for three kernels (stokeslet, stresslet, rotlet) of three-dimensional Stokes flow. The SE method computes the Fourier-space part of the periodic potential; it is based on the fast Fourier transform (FFT) and scales like O(N log N ) for N particles in the primary cell. Arbitrary periodicity (D = 3, 2, 1, 0 periodic directions) is supported, all within the same framework. This paper marks the completion and unification of the SE method for Stokes flow, uniting all kernels and periodic cases, thus allowing to perform efficient simulations of Stokes flow with arbitrary periodicity.
The modified kernels used to circumvent the singular behaviour of the original kernels form a crucial part of the method for D = 0, 1, 2. We have improved the convergence of the modified stokeslet and stresslet in the D = 0 case, and derived new modified kernels for the D = 1, 2 cases. With the new modified kernels, the truncation error is independent of periodicity, i.e. there is no penalty to the truncation error in the D = 0, 1, 2 cases compared to the D = 3 case. New improved truncation error estimates have been derived for the stokeslet and stresslet, valid in all periodic cases. We have also derived analytical formulas that can be used for validation for the D = 1, 2 cases.
Our version of the SE method uses the polynomial Kaiser-Bessel (PKB) window function introduced by [24], which reduces the support, and thus the computational time, needed to reach a given error compared with the truncated Gaussian (TG) window traditionally used in the SE method. We have showed numerically that the PKB window is indeed faster. The polynomial approximation of the window is selected such that it introduces no further error compared to the exact Kaiser-Bessel window. In the D = 1, 2 cases, an adaptive Fourier transform (AFT), first introduced by [52], is used to reduce the cost of the FFTs. In the D = 0 case, a precomputation step is used to reduce the computational time.
An automated procedure for selecting the parameters of the method has been presented and tested. The error is within one order of magnitude of the tolerance in the vast majority of cases. In the numerical results, we have furthermore showed that the pointwise error is well behaved, and that the computational time scales as expected. The method is fastest in the D = 3 case; the computational time for D = 2 is only slightly larger, while for D = 1 it is at most twice as large, and for D = 0 around three times as large as in the D = 3 case. The stresslet has more components than the other kernels and therefore requires around twice the computational time, due to more costly gridding and FFTs.
The SE method will be most efficient when applied to particles which are uniformly distributed in a primary cell with low aspect ratio (close to cubic). The computational time will increase if points are very unevenly distributed, or if the cell has a high aspect ratio. While the parameter selection procedure in principle works for any rectangular cuboid cell, it has been mainly tested for cells with low-to-moderate aspect ratios.
A drawback of the current parameter selection procedure is that only an absolute error tolerance, not a relative tolerance, can be given for the full potential (and real-space part). (The relative error can of course be computed a posteriori, but not a priori.) To be able to specify a relative tolerance for the full potential (and real-space part), one would need an estimate of its rms value. We have chosen not to construct such an estimate in this paper, since it would necessarily have to incorporate the distances from evaluation points to source points, making it more complicated than the Fourier-space part estimate that we constructed.
Previous iterations of the SE method have been successfully applied to simulations of Stokes flow based on boundary integral equations, in the D = 3, 0 cases. In the future, we will apply also the current version of the SE method to such simulations, with arbitrary periodicity. We believe that the techniques described in this paper would also be applicable to other kernels, such as the Rotne-Prager-Yamakawa tensor, assuming that a modified kernel can be derived (needed for reduced periodicity). Another possible future undertaking would be to examine the so-called "error pollution" and aliasing errors further and try to find a more refined rule, based on analytical formulas, to adjust the uniform grid spacing h and window size P .
The unified Spectral Ewald package with arbitrary periodicity includes the three kernels of Stokes flow as well as the harmonic kernel, and will be made available on GitHub [69]. that appears in (63), obtained by combining (62) and (61) for D = 2. Here, the kernel G is the stokeslet, rotlet or stresslet. As noted below (32)-(34), we can write G F = K A γ, where A is either the harmonic kernel (for the rotlet) or biharmonic kernel (for the stokeslet and stresslet). To be able to convert K to K, we introduce a new quantity (A.2) Inserting the relation G F = K A γ above, we find that This shows that in order to compute Q 2P , we can first compute the integral and then apply the relation To simplify the expressions in the following, we first introduce some auxiliary variables, namely (A.6) For the harmonic kernel H(r) = 1/|r| H(k) = 4π/|k| 2 and the Ewald screening function γ E (20), the integral Q H,2P , defined by (A.4), was stated by [60]. The value of the integral is found from [70, 3.954 (2), p. 504], with α and β + as in (A.6). Applying the relation (A.5) for the rotlet, with A = H and the operator K Ω given by (10), we find that the tensor Q Ω,2P is antisymmetric and given by Here, we have used that ∂β + /∂r 3 = αβ − . For the biharmonic kernel B(r) = |r| B(k) = −8π/|k| 4 and the Hasimoto screening function γ H (25), the integral Q B,2P is given by From the relation 10) between the integrands of (A.7) and (A.9), one obtains the relation Using the result (A.7) as well as ∂β + /∂α = r 3 β − − λ, we obtain with Λ as in (A.6).

Appendix A.2. Singular case
We now want to compute the integral that appears in (64), obtained from (62) and (61) for D = 2. Note that the integrand is now singular at κ 3 = 0. The integral does not exist in the Lebesgue sense, but we can interpret Q 2P,(0) as the (onedimensional) inverse Fourier transform of G F (0, 0, ·), in the distributional sense. (An introduction to the distributional Fourier transform can be found e.g. in [71, chapter 8].) In the same way as in Appendix A.1, we apply the relation G F = K A γ, noting that in this case, K(0, 0, κ 3 )e iκ3r3 = Ke iκ3r3 . The result is that where A is the harmonic kernel (for the rotlet) or biharmonic kernel (for the stokeslet and stresslet). Note that A γ = A F , so that One might try to compute Q A,2P,(0) (r 3 ; ξ) as the limit of Q A,2P (k 1 , k 2 , r 3 ; ξ) =:Ȃ F (k 1 , k 2 , r 3 ; ξ) from Appendix A.1 as α = k 2 1 + k 2 2 → 0. However, this limit exists neither for the harmonic nor for the biharmonic. Instead, we follow [60] and use the Ewald decomposition A = A R + A F to writȇ and this limit exists since the real-space part of the kernel has such a rapid decay in real space. It remains to computeȂ(k 1 , k 2 , r 3 ) andȂ(0, 0, r 3 ), but this is typically an easier problem since there is no screening function involved. For the harmonic, the derivation was done by [60]. The result is We now turn to the biharmonic kernel (A = B). Our first goal will be to computeB(0, 0, r 3 ) and B(k 1 , k 2 , r 3 ), which are needed in (A. 19)-(A.20). These are defined by taking the one-dimensional inverse Fourier transform in the free direction of B(k 1 , k 2 , κ 3 ) = −8π/(α 2 + κ 2 3 ) 2 , where again α = k 2 1 + k 2 2 . For α > 0, the inverse Fourier transform exists in the classical sense, and (using e.g. [72, 13.2, F41b, p. 319] and the convolution theorem) 3 ) 2B = 8πδ(r 3 ) for nonzero and zero α, respectively. We mention this since it corresponds to the view taken by [60] for the harmonic kernel.) We can now compute the limit in (A.20) to findB R (0, 0, r 3 ; ξ). WithB(k 1 , k 2 , r 3 ) taken from (A.23) andB F (k 1 , k 2 , r 3 ; ξ) taken from (A.12), we find that we want to compute the limit B R (0, 0, r 3 ; ξ) = lim Both the numerator and denominator have zero as their limits, so we apply L'Hôpital's rule. Differentiating the numerator with respect to α, we get (using that ∂Λ/∂α = −Λ/α − r 2 3 β + − αλ/(2ξ 2 )) πα −2r 2 3 e −α|r3| + r 2 3 β + + αλ 2ξ 2 , (A. 26) and the denominator becomes 3α 2 . Thus, the limit becomes This concludes the derivation for the biharmonic.
(In this case, i.e. D = 2, it is also possible to compute (A.30) directly by taking the one-dimensional inverse Fourier transform in the r 3 direction of B(0, 0, κ 3 ) γ H (0, 0, κ 3 ; ξ), in the distributional sense (e.g. using Wolfram Mathematica 12). However, the derivation presented above is more straightforward to apply in the D = 1 case. ) We now apply (A.16) for the stokeslet, with A = B and K S given by (9). The resulting tensor Q S,2P,(0) is symmetric and given by This result agrees with the expression for the k P = 0 mode derived by [20] using a different method. For the stresslet, with A = B and K T given by (11), we find that Q T,2P,(0) is given by for the inverse Fourier transform of a two-dimensional radial function f (κ 2 , κ 3 ) = f (κ). Here, ρ = r 2 2 + r 2 3 , κ = κ 2 2 + κ 2 3 , and J 0 (·) is the Bessel function of the first kind and order 0. Applying (B.27) to B yields B(k 1 , r 2 , r 3 ) = − 2ρ |k 1 | K 1 (|k 1 |ρ), k 1 = 0, (B.28) where K 1 (·) is the modified Bessel function of the second kind and order 1. Note that K ν (·) of one argument is not the same function as K ν (·, ·) of two arguments introduced in (B.8).
For k 1 = 0, the inverse Fourier transform must be interpreted in the distributional sense. By noting that the Fourier transform of − (∇ 2 2 + ∇ where B is an arbitrary positive gauge constant, and we will here set c B = 0. We now compute the limit found in (B.22), withB(k 1 , r 2 , r 3 ) taken from (B.28) andB F (k 1 , r 2 , r 3 ; ξ) taken from (B.14). The limit is Recall that k 1 is found also in U , see (B.6). Using the relation [73] K ν (a, b) + K −ν (b, a) = 2(a/b) ν/2 K ν (2 At this point, we note that −ρ 2 log(ρ/ B ) = −ρ 2 log(ρ) + ρ 2 log( B ) and that the contribution to Q S,1P,(0) from the term ρ 2 log( B ), when applying K S for the stokeslet, given by (9) where K 3 is given by (40), and S F is given by (35). We will assume that (C.1) is truncated outside a sphere of radius |k| = k ∞ for some k ∞ > 0. (In the SE method, the sum is actually truncated outside a cube of side length 2k ∞ , but that can only lead to a smaller error, and the difference will in any case not be large.) The pointwise error from truncating (C.1) is given by To be able to follow [23], the sum over k is approximated by an integral. Multiplying and dividing (C.2) by (∆k) 3 , where ∆k = 2π/L is the wavenumber resolution, allows us to make the approximation e S,F (x; ξ) ≈ 1 (2π) 3 such that e S,F (x; ξ) ≈ N n=1 E (x − x n )f (x n ). As done by [23], we approximate the directional component by its root mean square value (computed using spherical coordinates), and keep only the highest-order term with respect to |k| in (C.4), since that term will dominate the error for large k ∞ . Introducing spherical coordinates [k, θ, ϕ], with the coordinate system chosen such that k · r = kr cos θ (where k = |k| and r = |r|), we get the result E jl (r) ≈ 2π √ 2 12π 2 ξ 2 ∞ k∞ π 0 e −k 2 /(2ξ) 2 e ikr cos θ k 2 sin θ dθ dk (C. 6) where the integral over ϕ could be evaluated to 2π directly. Computing also the integral over θ yields The exponential decay of e −k 2 /(4ξ 2 ) makes sure that the dominant contribution comes from the beginning of the interval, where k ≈ k ∞ . This allows the approximation E jl (r) ≈ √ 2 3πξ 2 k ∞ r ∞ k∞ sin(kr)e −k 2 /(2ξ) 2 dk. (C.8) The remaining integral can be computed exactly in terms of the error function, but to get a more manageable error estimate, we will approximate it. Writing sin(kr) = Im{e ikr }, and once again using k ≈ k ∞ , we get ∞ k∞ sin(kr)e −k 2 /(2ξ) 2 dk = Im ∞ k∞ e ikr−k 2 /(2ξ) 2 dk = Im ∞ k∞ ir − k/(2ξ 2 ) ir − k/(2ξ 2 ) e ikr−k 2 /(2ξ) 2 dk ≈ Im 1 ir − k ∞ /(2ξ 2 ) ∞ k∞ ir − k/(2ξ 2 ) e ikr−k 2 /(2ξ) 2 dk = e −k 2 ∞ /(2ξ) 2 r 2 + k 2 ∞ /(4ξ 4 ) r cos(k ∞ r) + k ∞ 2ξ 2 sin(k ∞ r) . (C.9) Again assuming that k ∞ is large (also compared to r), we keep only the highest-order terms with respect to k ∞ . Inserting (C.9) into (C.8) then yields sin(k ∞ r) r e −k 2 ∞ /(2ξ) 2 =: E(r), where V is a volume containing all vectors r mn = x m − x n , and Q is defined as in (154). Selecting V to be a sphere of radius ρ := √ 3L/2, we get (inserting (C.10) and computing the integral over V using spherical coordinates) Since k ∞ is large, the term involving sin(2k ∞ ρ)/k ∞ can be neglected. This leads to the estimate which is (155).

Appendix C.2. Stresslet estimate
The derivation for the stresslet (37) proceeds in exactly the same way as for the stokeslet, with the only difference being that the stresslet has an extra factor of |k| compared to the stokeslet, and that the root mean square of the directional component is instead of √ 2/3. This means that the stresslet estimate will differ from the stokeslet estimate by a factor of k ∞ 7/6. Multiplying (155) by this factor yields (156).