{'fraction_non_alphanumeric': 0.03357366665543614, 'fraction_numerical': 0.010410700448098109, 'mean_word_length': 4.229759492555722, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A fast integral equation method for solid particles in viscous flow using quadrature by expansion

Boundary integral methods are advantageous when simulating viscous flow around rigid particles, due to the reduction in number of unknowns and straightforward handling of the geometry. In this work we present a fast and accurate framework for simulating spheroids in periodic Stokes flow, which is based on the completed double layer boundary integral formulation. The framework implements a new method known as quadrature by expansion (QBX), which uses surrogate local expansions of the layer potential to evaluate it to very high accuracy both on and off the particle surfaces. This quadrature method is accelerated through a newly developed precomputation scheme. The long range interactions are computed using the spectral Ewald (SE) fast summation method, which after integration with QBX allows the resulting system to be solved in M log M time, where M is the number of particles. This framework is suitable for simulations of large particle systems, and can be used for studying e.g. porous media models.


Introduction
Fluid flows involving microscopic, rigid particles are common both in nature and in industrial processes. The macroscopic properties of such systems are often determined by the particle interactions happening on the smallest scale of the flow. Examples of such systems can be found in sedimenting suspensions [17] and porous media modeling [13]. To fully understand the interactions between particles, and the effect those interactions have on the fluid flow, numerical simulation is a valuable tool. Stokes equations are often valid in this context, due to the small particle size and low fluid velocity. For problems governed by Stokes equations it is possible to use boundary integral methods, where the solution is represented as a layer potential from the boundaries of the domain (in this case the particle surfaces). Since the problem then is formulated as a boundary integral equation on the union of the particle surfaces, the dimension of the domain which has to be discretized is reduced from R 3 to R 2 , significantly reducing the number of unknowns. This also removes the problem of combining a volume grid with moving boundaries, which can be challenging. However, boundary integral methods come with a set of challenges of their own, two of which can be considered major.
The first major challenge of boundary integral methods is that the resulting linear system after e.g. a Nyström discretization is dense, such that the cost of one left hand side evaluation is O(N 2 ), where N is the number of unknowns. This puts a severe limit on the size of the systems that can be considered, even if a rapidly converging interative method is used. This can be overcome by evaluating the layer potential using a fast method, such as the fast multipole method (FMM) [16] or a fast Ewald summation method [14,23] (for periodic problems). These methods reduce the cost of one left hand side evaluation to O(N ) and O(N log N ), respectively, thereby making boundary integral methods suitable for large-scale computations.
The second major challenge of boundary integral methods is that one needs accurate quadrature of singular and nearly singular integrals when evaluating the layer potential. Developing such methods which are accurate, fast and work for arbitrary geometries is a topic of current research, particularly for the problem of nearly singular quadrature in three dimensions. In two dimensions there are efficient methods for nearly singular quadrature which gain a lot of their power from the complex variable formulation [18,6]. In three dimensions the situation appears less resolved, though several different methods have been successfully used in practical applications [9,10,36,28,33,37]. A common feature of many of these methods is however that they are highly target specific, meaning that the cost grows rapidly if there are many nearly singular integrals to be evaluated.
Quadrature by expansion (QBX) [20,7] is a fairly new method for numerical integration of singular and nearly singular integrals. The method is built around evaluating layer potentials through local expansions, and comes equipped with a solid convergence theory [15]. It was originally presented for the Helmholtz kernel in two dimensions, but the principles of the method generalize directly both to three dimensions and other kernels. A promising feature of QBX is that it is possible to combine it with an FMM, thereby creating a fully O(N ) method that is able to accurately evaluate the potential everywhere.
In this paper we extend the QBX framework to deal with the Stokes double layer potential in three dimensions, and combine it with the fast Ewald summation method presented in [2]. The result is a robust and scalable framework for computing Stokes flow around periodic systems of rigid, spheroidal particles. To limit the scope of the present paper, we restrict our attention to stationary particles. The method could easily be extended to dynamic problems -such as sedimentation -by coupling it to an ODE solver, as was done in [2].
The structure of this paper is as follows: In section 2 we state the necessary boundary integral formulation for Stokes flow around rigid particles. In section 3 we introduce QBX and discuss errors and parameter selection. In section 4 we show how QBX can be accelerated for our problem, leading to a computationally feasible method. In section 5 we briefly touch on the subject of Ewald summation, and how to combine it with QBX.
Finally, in section 6 we present selected numerical results, which we draw both from validation tests and from an example application. We also include an appendix (A), which covers the computational details of the acceleration scheme presented in section 4.

Formulation
Our application of interest is that of rigid particles in a Newtonian fluid. The particles are assumed small enough that the Reynolds number is close to zero, such that we can model the flow as being Stokes flow, governed by the Stokes equations, We here limit ourselves to spheroidal particles with semi-axes a and c. The surface of such a particle, oriented along the z coordinate axis, can in Cartesian coordinates be described as If parametrized in the spherical coordinates θ ∈ [0, π] and ϕ ∈ [0, 2π), the same surface is given by x = a sin(θ) cos(ϕ), y = a sin(θ) sin(ϕ), z = c cos(θ). ( The particles can be classified into three distinct subgroups: prolate (a < c), spherical (a = c), and oblate (a > c).

Boundary integral formulation
We are considering Stokes flow as described by the Stokes equations (1), which by virtue of being a set of linear partial differential equations with constant coefficients have solutions that can be represented using boundary integrals. For these representations, the essential Green's functions are the stokeslet S, the stresslet T and the rotlet R 1 .
where r = x − y and r = |r|.
To compute the flow associated with a rigid particle in a viscous fluid, we use the completed boundary integral representation of Power & Miranda [26], where the flow in the domain is expressed as a double layer potential D plus a completion flow V and a background flow u bg , The double layer potential is the flow associated with a double layer density q, defined on the particle surface Γ. It is computed by integrating q over Γ together with the stresslet singularity T and the outward unit normal n, The integral in (8) is weakly singular for x ∈ Γ, and is then taken as the principal value integral, assuming that Γ is Lyapunov smooth. Letting x approach the surface from either the interior or exterior domain, there is a jump in the double layer potential, The completion flow V is added to the formulation due to the double layer potential's inability to represent a net force f and torque t on the particle [26]. For a rigid particle, a suitable completion flow is that generated by a stokeslet singularity S of magnitude f /8πµ and a rotlet singularity R of magnitude t/8πµ, both located at the particle center x c , For a particle undergoing rigid body motion with translational velocity U and rotational velocity Ω, we set a no-slip boundary condition, Letting x in our flow field representation (7) go to the surface of the particle from the exterior, we get a diagonal term from (9). We then enforce the no-slip boundary condition (11) on the surface, and close the system using (16) and (17). This results in a boundary integral equation (BIE) of the second kind in the density q, The fact that (12) is second kind with a compact integral operator is the main benefit of using a double layer formulation. After discretization using a quadrature rule and the Nyström method, the corresponding linear system can be said to be well conditioned in two senses; one is that the condition number is bounded (and usually very low), the other is that the condition number stays constant under grid refinement (this is in contrast to single layer formulations, where the condition number increases under refinement). These properties make the system suitable for iterative solution using the generalized minimum residual method [30] (GMRES), which typically converges rapidly.
The equation system in (12) needs to be closed by specifying either the rigid body motion (U , Ω) or the completion flow V. This corresponds to using two different formulations: the resistance problem formulation or the mobility problem formulation.

The resistance problem
When the rigid body motion (U , Ω) of the particle is known, the problem is to find the resulting force f and torque t exerted by the fluid on the particle. These can be related to the double layer density q through the formulation by Power & Miranda [26], such that the completion flow is a functional of q,

The mobility problem
When the external forcing (f , t) on the particle is known, the problem is to find the rigid body motion (U , Ω) of the particle. Using the formulation available in e.g. Pozrikidis [27], the rigid body motion vectors can be computed as functionals of q, where S Γ is the surface area of Γ and The vectors ω (n) , n = 1, 2, 3, are independent unit vectors satisfying

Particle systems
For a system of M particles, the linearity of the Stokes equations allows us to express the velocity field as a superposition of the fields from the individual particles, where Γ α is the surface of particle α, and (V α , U α , Ω α ) are defined either as in the resistance formulation or the mobility formulation.
Evaluating (20) and enforcing the no-slip boundary condition (11) on the surface of each particle gives the boundary integral equation for a system of particles,

Discrete formulation
To solve the (single particle) BIE (12), we first need a way of numerically evaluating integrals over the particle surface, which we denote This we do using a quadrature rule Q N , which defines a set of N nodes x i and weights w i on Γ, such that The details of the quadrature used will be discussed further in section 3. We denote by the superscript h a quantity computed using Q N , e.g.
Applying the Nyström method, where the integral equation is enforced at the quadrature nodes, the BIE for the single-particle system (12) is approximated by the 3N × 3N linear system in the density values q(x i ), where either V = V h or U = U h and Ω = Ω h , depending on if the problem is a resistance problem or a mobility problem.
The matrix corresponding to the linear system (26) is dense, so the cost of applying a direct solution method to it would be O(N 3 ). However, when solving with GMRES the number of iterations required to reach a certain tolerance is largely independent of the discretization. The cost is then O(N 2 ), with a constant that depends on the problem geometry. This is acceptable for a single particle, since the solution q usually can be resolved with a limited number of surface nodes. For a problem with multiple particles, however, the corresponding linear system (derived from (22)) is 3M N × 3M N and dense. The cost then grows quadratically with the number of particles, which severely limits the size of the problem that can be considered. The remedy for this is to use a fast method, which exploits the structure of the problem to evaluate the left hand side of the problem (also known as the matrix-vector product, or matvec) in less than quadratic time. Common fast methods are the FMM and, for periodic problems, fast Ewald summation. The respective complexities for these two methods are O(M ) and O(M log M ) as the number of particles is increased. We have here used fast Ewald summation, which will be further discussed in section 5.

Quadrature
To derive the linear system (26) corresponding to the BIE (12), we needed a quadrature method Q N (24) for evaluating the double layer potential D (8). We will in this section describe such a method, which is able to accurately resolve D also when the integral is nearly singular or singular.

Spheroidal grid
When performing quadrature on a spheroidal particle, it is natural to use the parametrized description (3) in (θ, ϕ) and to discretize it using an n θ × n ϕ grid in those coordinates, with the grid aspect ratio given by This gives us a total of N = n θ n ϕ quadrature points on the surface. We let (θ i , λ θ i ), i ∈ {1 . . . n θ }, be the nodes and quadrature weights of an n θ -point Gauss-Legendre quadrature rule on the interval [0, π], and (ϕ j , λ ϕ j ), j ∈ {1 . . . n ϕ }, be the nodes and quadrature weights of the trapezoidal rule on the interval [0, 2π). This gives us the quadrature rule (a) Stokes flow around a single translating spheroid, exact solution from [12].
(b) Error when evaluating double layer potential from spheroid using direct quadrature. Figure 1: Example of error due to nearly singular integration (b), when evaluating flow due to a single particle (a).
where W (θ, ϕ) is the area element. For ease of notation we can simplify this expression by letting one index cover all quadrature points and introducing y i and w i such that We denote this grid and the accompanying quadrature rule the spheroidal grid, and will refer to the quadrature rule as the direct quadrature. The direct quadrature rule defined in (28) allows us to approximate integrals over the surface with spectral accuracy for smooth and well-resolved integrands, and works well for evaluating the double layer potential (8) when the target point x is far enough away from the surface Γ. However, when x is on Γ the integral is singular and can not be computed using a quadrature rule for smooth functions. There is also a difficulty when x lies in the domain, but is close to Γ. The singularity in the stresslet then causes the integral kernel to be sharply peaked, and trying to evaluate it using a smooth quadrature rule gives an error which grows exponentially as x approaches Γ (see Figure 1). We refer to this as the integral being nearly singular.
The singular case (x ∈ Γ) can for the double layer potential be evaluated using a singularity subtraction method, which gives third order accuracy at virtually no increased cost [2]. For higher order accuracy in the singular case, and for evaluating the nearly singular case, specialized quadrature methods are required.
In what follows of this section we will first give a brief introduction to an existing quadrature method for singular and nearly singular layer potentials that is based on local expansions. We will then show how this method can be used to evaluate the Stokes double layer potential. Lastly, we will show how this method in our case can be accelerated by geometric considerations.

Quadrature by expansion
Quadrature by expansion (QBX) is a new quadrature scheme first introduced for the Helmholtz equation in two dimensions by Klöckner et al. [20]. The method is based on the observation that a straightforward quadrature rule applied to a layer potential is inaccurate when evaluated close to or on the surface, but well resolved when sufficiently far away, as illustrated in Figure 1. Since the layer potential is a smooth function away from the boundary, we can pick any point c in the domain and create a local expansion of the potential about that point. Picking that point in the well-resolved part of the domain, we can compute the coefficients of the expansion using our direct quadrature rule. The domain of convergence for the local expansion is the ball of radius r, centered on c, and so the expansion can be used to accurately evaluate the potential close to the surface and at the point where the ball touches the surface (see Figure 2 for an example geometry). This is the essential idea of QBX. For further reading regarding the convergence and evaluation of these local expansions, we refer to Epstein et al. [15] and Barnett [7]. In three dimensions, a good starting point for the construction of local expansions for layer potentials is the expansion of the Green's function for the Laplace equation, called the Laplace expansion, where Y m l is the spherical harmonics function of degree l and order m, and (r x , θ x , ϕ x ) and (r y , θ y , ϕ y ) are the coordinates of x and y in a spherical coordinate system centered at c. This expansion gives us the separation of source (y) and target (x) that is necessary for the method. As a first application, we will consider the double layer potential of the Laplace equation, which we henceforth will refer to as the dipole potential. For a smooth vector density ρ defined on the surface Γ, we define it as Inserting the expansion (31) into (32) and moving the terms related to the target point x out of the integration, we get where z lm are the local expansion coefficients of the potential, Once the expansion coefficients z lm are known, they can be used to evaluate the potential at any target point x in the domain of convergence. At this stage we wish to highlight two things: (i) We have so far made no approximations, evaluating the potential through expressions (33) and (34) is still exact. The sum in l will need to be truncated, and we will discuss the truncation error that it introduces in section 3.5. (ii) The integrand in (34) contains no singularity, and can be evaluated using a smooth quadrature rule. The quadrature errors related to this smooth quadrature will also be discussed in section 3.5.

QBX for the stresslet
To evaluate the Stokes double layer potential (8) using QBX, we need to be able to expand its kernel in a way that separates source and target, just as we could for the dipole potential. One way of doing this is to express the double layer potential in terms of the dipole potential, for which we already have the expansion (33). Here we use the result by Tornberg and Greengard [34], which was developed with the purpose of using a harmonic FMM to compute the flow due to a collection of stokeslets and stresslets. Denoting by F the kernel of the dipole potential, the kernel of the double layer potential can be expressed as Rearranging terms, the integrand of (8) can then be written Consequently, the double layer potential can be expressed in terms of dipole potentials as With this result we can create an expansion of the double layer potential by creating four expansions of the dipole potential with four different densities (including the implicit summation over j). This means that the three vector components of the double layer potential are computed using four scalar expansions.The expansion coefficients that need to be computed for a given expansion center c are then The spherical harmonics are conjugate symmetric in m, Y −m l (θ, ϕ) = Y m l (θ, ϕ) * , and since the density q is real the same conjugate symmetry holds for the expansion coefficients, Hence, we only need to compute the coefficients for 0 ≤ m ≤ l. Truncating the expansion at some p = l max , the total number of coefficients that we need to compute for each component j is

Numerical scheme
To solve the discrete BIE for a single particle (26), we need to evaluate the double layer potential at the quadrature points x 1 , . . . , x N . A local expansion of the potential is only convergent at a single point on the surface, where the ball of convergence touches the surface, so we need N expansion centers c 1 , . . . , c N located at a distance r normally from the quadrature points, The same expansions are used to evaluate the potential on the surface and in the domain. A natural choice of expansion radius is then where d has been determined as the minimum distance from the surface at which the direct quadrature can be used for a given error tolerance . For target points x such that min y∈Γ |x − y| < d , the potential is then evaluated using the closest expansion center.

Error analysis
Here we give a sketch of an error analysis for the dipole potential, with the argument that it carries over to the double layer potential, as the latter is constructed as a combination of the former. For more in-depth discussions we refer to [15] and [3]. When evaluating the dipole potential (32) using QBX, we truncate the expansion (33) at some order p and use our direct Gauss-Legendre/trapezoidal quadrature rule (28) to approximate the coefficients z lm (34). Denoting the approximate coefficientsz lm , the error is Adding and subtracting the exact coefficients z lm to the expansion, the error can be split into two components, where E T is the truncation error due to using a finite number of terms in the expansion, and E Q is the quadrature error due using a discrete quadrature rule, From the results in [15] we have that the truncation error can be expected to behave as In [3] an estimate was derived for the QBX quadrature error, when considering the Laplace single layer potential on a spheroid. Simplifying those results, we can expect that the quadrature error should behave approximately as for some α, β > 0 and h a measure of the grid size; for the spheroidal grid we define Putting the two expression together, we see that the total error has one component that that decays with p and one component that grows with p, This means that if we fix r and h for a given problem, then there exists an optimal value for p that minimizes the error. If we instead fix the relation r/h at some constant value, e.g. r = 3h, then the quadrature error stays constant (E Q = C Q ) such that for a given p The scheme can then be viewed as having order of accuracy (p+1) under grid refinement, up to the point where the truncation and quadrature errors are equal. After that we have to refine the grid (while maintaining a constant r) to get any additional improvements.
In practice we want to be able to get arbitrarily high accuracy out of QBX, without having to introduce more degrees of freedom. To overcome the problem of the quadrature error growing with p, we therefore upsample the surface density to a factor κ finer grid before computing the coefficients. The fine grid is of the same type as the original, but with κn θ × κn ϕ points, and the density on the fine grid is computed using trigonometric interpolation and barycentric Lagrange interpolation [8]. This allows us to write the total error as The truncation and quadrature errors can then be controlled separately, since the exponential term will rapidly reduce the quadrature error when κ is increased.
To illustrate the effects of the parameters p and κ, we solve the mobility problem for a single oblate spheroid (a/c = 2) with a 20 × 40 grid using r/h = 3/2 and a range of p and κ. Measuring the L 2 norm of the error in the density compared to a reference solution, we can see in Figure 3a that the error decreases exponentially in p until it hits a floor given by the resolution of the underlying grid, but only for κ large enough. If κ is taken too small, the error increases almost exponentially after reaching a minimum. Finding the optimal κ for a given p is not trivial, and picking κ larger than necessary is costly, since it adds a factor κ 2 to the computational complexity.

Two-sided expansions
When evaluating the singular double layer potential on a surface, we want the principal value integral. However, what we get when evaluating it using QBX is the one-sided limit from the side of the expansion center, which we denote D + or D − depending on side,  The relation between the principal value and the one-sided limit is known through the jump relation (9), which can be used together with the QBX result to recover the principal value. However, and this was originally noted in [20], it turns out that the spectrum of the linear system better approximates that of the original operator if we compute the principal value through the relation To do this in practice, we need to put one expansion center on each side of the target point on the surface and take the average potential from the two to get the principal value, as illustrated in Figure 3b. This seemingly doubles the amount of work needed, but through the acceleration scheme of section 4 that extra cost is hidden from the computations.

Parameter selection
For a given geometry, grid resolution h and tolerance , an empirical process for selecting the QBX parameters p, κ and r is as follows: 1. Set r = d /2, such that the error from the regular quadrature will be less than epsilon when the evaluation point is at least a distance d from Γ.
2. Set p large enough such that the truncation error E T from the QBX expansion is smaller than for points closer than d from Γ.
3. Set κ large enough to match p, such that the quadrature error E Q is kept below . Figure 4 shows an example where the QBX parameters have been tuned for high accuracy.   , 2), compared to exact solution from [12]. Grid is 60 × 30, QBX parameters are p = 40, k = 10, r/h = 2.
Note the error pattern inside the QBX area, which is due to an error in the solution, not the quadrature.

Acceleration of QBX on spheroidal grid
Evaluating the potential from one particle using the QBX scheme described in section 3 involves a lot of work. For a spheroidal grid with N points, a given upsampling factor κ and expansion order p, we need to interpolate the density to a factor κ 2 larger grid and then compute O(p 2 ) expansion coefficients at every expansion center. This puts the cost for computing the potential from a particle on itself, the "self-interaction", at O(κ 2 p 2 N 2 ). One natural way of speeding up the process is to incorporate the QBX scheme into an FMM, which hierarchically computes potentials using expansions of the same kind as those used in QBX (the method was in fact conceived with this in mind). This makes it possible to compute everything on the fly at a computational cost that is linear in N (but still includes κ 2 p 2 ). The modifications are however nontrivial, and the only unified FMM/QBX method published to date (by Rachh et al. [29]) is for two-dimensional problems.
Another way of achieving speedup is to precompute as much as possible of the work associated with QBX, and then use precomputed values whenever possible. This comes at a potentially large storage cost, and still has a computational complexity that is quadratic in N , albeit with a significantly smaller leading constant. In this work we have used precomputations in combination with a fast Ewald summation method, which computes the long-range interactions between M particles in O(M log M ) time (more on this in section 5). This is feasible partly because the number N of grid points required on a spheroid is limited, such that the overall computational cost scales like O(M log M ), and partly because the spheroidal geometry of the particles in our problem significantly reduces the storage required for the precomputed values. We will in this section outline the principles and main results of our precomputation strategy; the details of the implementation are available for reference in appendix A.

Self-interaction
When solving the boundary integral equation (12) on a single particle using the spheroidal grid of section 3.1, we typically use the Nyström method to create a linear system which we solve iteratively. We then have N grid points denoted x 1 , . . . , x N on the surface, and for every iteration we need to compute the self-interaction at those points. With The matrix R i then represents the resulting action after upsampling the density to a fine grid, using the fine grid to compute the expansion coefficients at the expansion center c i , and finally evaluating that expansion at x i . Precomputing all the matrices R 1 , . . . , R N requires 3N × 3N memory storage, but once computed they allow the computation of Q → u(x i ) at all points x i at a the cost of a matrix-vector product, i.e. O(N 2 ). This means that the O(p 2 κ 2 ) cost related to the QBX parameters -and ultimately to the accuracy -only appears in the precomputation step. Once the R i matrices are computed the self-interaction evaluation has a fixed cost, independent of which QBX parameters were used. The memory cost of precomputing all the self interaction matrices can be reduced by taking advantage of the fact that the spheroidal grid is rotationally symmetric about the polar axis. This symmetry means that the matrices R a and R b related to two points x a and x b on the same latitude (θ a = θ b ) can be related as where T z ∈ R 3×3 represents a rotation about the polar axis of the particle, and P ∈ R N ×N is a permutation matrix, i.e. a row permutation of the identity matrix that represents a permutation of the point ordering. The symbol ⊗ denotes the Kronecker product. As a consequence, it is enough to precompute matrices related to the n θ points on the first longitude, and then use rotations and permutations to compute the potential at all grid points with those matrices. The savings in memory and precomputation time can be further increased by taking advantage of the fact the spheroidal grid has a mirror symmetry about the equator. For two points x b and x c on the same longitude (ϕ b = ϕ c ) but different sides of the equator, such that θ b = π − θ c , we can relate the matrices R b and R c as where T x ∈ R 3×3 represents a rotation that mirrors the points x b and x c into each other, and F ∈ R N ×N is a permutation matrix. Taking advantage of this mirror symmetry (a) Spheroidal grid with the n θ /2 first points circled.
(b) Example of configuration where nearsingular quadrature is necessary. Figure 5 allows a further reduction of the required number of precomputed matrices to n θ /2, corresponding to the points marked in Figure 5a.

Nearly singular evaluation
When particles are close to each other, such as in Figure 5b, we need to use QBX to evaluate their interaction accurately, since at least some of the quadrature points on one particle will be too close to the other particle for the direct quadrature to be accurate.
To be able to evaluate the potential at arbitrary points close to a particle using QBX, we need the expansion coefficients for the potential at all the expansion centers in the domain outside the surface. We let z j (c i ) ∈ C Np , j ∈ {1 . . . 4}, denote the four vectors of N p (41) coefficients each at a an expansion center c i , as defined in (39). Analogously to (56), the process of upsampling the density and computing the expansion coefficients at all the expansion coefficients can then by represented by a set of matrices M j i ∈ C Np×3N , such that Using the same arguments of symmetry as above, we can relate these matrices by rotational symmetry as and by mirror symmetry as where {Θ, E} ∈ R Np×Np are diagonal, and {A, B} ∈ R 4×4 orthogonal. Just as for the matrices R i , it is therefore sufficient to precompute the matrices M j i for the first n θ /2 expansions points, since the remaining matrices can be reconstructed by symmetry operations.
The precomputation of the matrices M j i hides the cost κ 2 related to the upsampling from the nearly singular evaluation, though the O(p 2 ) expansion order cost still affects the evaluation and storage cost. This is because it reflects the number of expansion coefficients used, and we need all the coefficients to able to evaluate the expansion at arbitrary target points. Only when precomputing for a specific target point, as in the self-interaction, can we get rid of the p 2 factor.

Storage complexity
Precomputing the density-to-coefficient matrices M j i for a single expansion center c i requires storage of 12N p N complex numbers, while the density-to-potential matrix R i related to a target point x i only requires storage of 9N real numbers. Considering a spheroidal body, where N = n ϕ n θ , and taking both symmetries into account, the storage required for computing all necessary M j i is then 6N p n θ 2 n ϕ complex numbers, while the storage required for all necessary R i is 9 2 n θ 2 n ϕ real numbers. As an example, using a set of 50 × 50 quadrature points and p = 15 requires approximately 1.6 GB for the M j i and 4.5 MB for the R i , using double precision. The number of expansion coefficients scales as N p = O(p 2 ), so assuming n ϕ ∝ n θ = n allows us to write the scaling of these expressions as O(p 2 n 3 ) and O(n 3 ). Alternatively, introducing a grid size h ∝ 1/n gives O(p 2 h −3 ) and O(h −3 ).

Summary of precomputation scheme
The precomputation scheme outlined above minimizes the cost for QBX when computing it directly, as opposed to accelerating it by an FMM. For a given setup of particle shape, quadrature grid and QBX parameters, we can precompute the matrices R i and M j i , associated with self-interaction and nearly singular evaluation, and store them to disk. The symmetries of the spheroidal grid allow us to do this at a manageable cost in terms of memory. Once we have these precomputed values, we can load and use them for simulations with a large number of particles that all have identical shape and grid. The single set of precomputed values can then be used for all particles in the simulation, since they will differ from each other only by a rigid body rotation.
For the self-interaction, the costs of upsampling and expansion order are hidden from the evaluation step by the precomputations. We can therefore pick κ and p large enough to make sure that that the accuracy in the evlauation of the double layer potential is only limited by how well the density q is resolved on the underlying quadrature grid. For the nearly singular evaluation the expansion order p still affects the computational cost, so there is still a cost/accuracy tradeoff there. However, the upsampling factor κ is hidden by the precomputations, so we can always make sure that we pick it large enough to resolve integrands associated with large p.

Periodic systems and Ewald summation
To be able to model a large particle system without boundary effects, we use the common method of applying periodic boundary conditions. We then consider a system of M particles contained in a box of size L 1 × L 2 × L 3 , which is periodically replicated in all three spatial directions. The velocity field is then periodic, where the superscript P denotes the periodicity and τ is a periodic shift, The periodic velocity field is computed as a superposition of the field from all periodic replications of all particles, where D P is the periodic double layer potential, and V P,α is the periodic completion flow, Computing u P by direct summation over all periodic images is not feasible due to the sums over τ in (65) and (66) only being conditionally convergent. Here we compute them using a fast Ewald summation method, which we have also used in previous work [2]. The stresslet T in (65) is then decomposed into two parts, where T (R) is sharply peaked and short-ranged, while T (F ) is smooth and long-ranged. The Ewald parameter ξ determines how short-ranged and smooth the two respective parts are. Away from the source y the short-ranged part T (R) decays exponentially, This makes it possible to introduce a truncation radius r c and only consider near-neighbor interactions within that radius (|x − y| ≤ r c ). The long-ranged part T (F ) converges rapidly when computed in Fourier space, where its Fourier coefficients T (F ) (k, ξ) decay exponentially due to its smoothness, The Fourier expansion is truncated at a maximum wave number K = max k k ∞ and computed using the spectral Ewald (SE) method [22,23], which uses a fast Fourier transform (FFT) to compute the long-range interactions to spectral accuracy. In [2] we describe in detail how this is done for the stresslet, and provide accurate truncation error estimates which can be used for selecting the Ewald parameter ξ and the cutoffs r c and K. The equivalent development for the stokeslet and rotlet singularities are available in [22] and [1]. Computing T (F ) using the SE method and only computing T (R) for near neighbors gives a method with O(N log N ) complexity, where N is the total number of degrees of freedom when the number of particles in the system is increased, under the assumption that the average number of near neighbors of each discrete point is kept constant [2].

QBX and Ewald summation
In section 3 we developed a framework for computing the double layer potential (8) to high accuracy for both the singular and nearly singular cases, based on an expansion for the stresslet T = r i r j r k /r 5 . When computing the periodic double layer potential (65) using Ewald summation, we have split the stresslet into a short-ranged part T (R) and a long-ranged part T (F ) (67). The explicit form of the short-ranged part is [2] T (R) where r = x − y and r = |r|. All of the singular behavior of the stresslet is contained in T (R) , since a series expansion of T (R) in the limit of r → 0 shows that This means that we need a way of using QBX for evaluating the short-ranged potential, which we define as ijk (x, y, ξ)q j (y)n k (y)dS y .
Deriving a local expansion for T (R) , as defined in (70), is a daunting proposition. Instead, we rewrite the short-ranged potential as The first term is now the ordinary double layer potential (8), for which we can use our existing QBX framework. The second term is the smooth and long-ranged part of the potential, which we can evaluate using our direct quadrature rule. In our implementation we use (74) together with QBX if x is close to or on Γ, and (72) truncated at r c otherwise.

Computational complexity
When now consider a system of M particles with N degrees of freedom each, which we solve using our accelerated QBX scheme combined with SE. The computational complexity of solving this system then has three distinct components: (i) the interactions of particles with themselves (the self-interactions), (ii) the nearly singular interactions between particles, and (iii) the well-separated interactions. The self-interactions are completely precomputed, and can therefore be evaluated at a O(M N 2 ) cost. The nearly singular interactions use the precomputed M j i -matrices, and contribute with a O(p 2 M N 2 ) cost, p being the QBX expansion order and the leading constant depending strongly on the interparticle spacing. The well-separated interactions are handled by SE, at a cost which scales as O(M N log M N ) as we increase the size of the particle system. This puts the total computational complexity at (75) For a given particle system with fixed N , p and average interparticle spacing, the cost then grows as O(M log M ) as the system size is increased by increasing the number of particles M .

Numerical results
We can now summarize our numerical scheme as follows: • Spheroidal grid: For a spheroidal body with axes a and c, we introduce a spheroidal grid (section 3.1) of dimensions N = n θ × n ϕ . For a system with M particles, an identical grid is used for each particle.
• QBX: We place expansion centers c i (42) at a normal distance r from the grid points x i , and use an order p local expansion (section 3.2) to evaluate the double layer potential (8) in the vicinity of each expansion center. The expansion is computed using a factor κ upsampled grid. The QBX computations are accelerated by the precomputation scheme of section 4.
• Ewald: For periodic problems, we compute the velocity field by combining QBX with a fast and spectrally accurate Ewald summation method (section 5).
• Solution: We solve the discrete system (26) with 3M N unknowns iteratively using GMRES, with the left hand side computed using QBX and Ewald.
• Preconditioning: For systems with many particles, we construct a simple blockdiagonal preconditioner by computing the explicit inverse of the corresponding single-particle system and using it to create the diagonal blocks (through rotation). This requires (3N ) 2 memory storage, but reduces the number of iterations by as much as a factor three for some problems.
This scheme allows us to efficiently compute the Stokes flow around systems of particles where the particles are very close to each other, and also to accurately evaluate the flow everywhere in the domain. An example is shown in Figure 6, where one can clearly see the benefit using QBX when evaluating the potential close to the particles.

Validation
To validate that our method is correct, we have compared it against several analytical, semi-analytical and numerical results. For a single particle in free space, analytic results are available in Chwang & Wu [12,11] and Jeffery [19]. Figure 4 shows an example of a comparison against [12]. For systems of three spheres in free space, semi-analytic results for several configurations are available in Wilson [35]. For a periodic array of spheres, numeric results are available in Zick & Homsy [38] and Sangani & Acrivos [31]. The described method can solve all these problems to arbitrary accuracy, provided that the grid is fine enough to represent the solution, and that the QBX parameter selection of section 3.7 is followed. We will below discuss results for a three-sphere system and for a periodic array of spheres. Our implementation of the method has also been used in a separate work by Bagge [5], where it was thoroughly validated against the solutions by Chwang [11] and Jeffery [19].

Triangle of spheres
To test the ability of our method to resolve particle interactions in the lubrication limit, we set up the following test case from Wilson [35]: Three spheres of radius a are arranged in an equilateral triangle with sides sa, s > 2. A force of magnitude 6µπa is applied to sphere 1 towards the centroid of triangle, pushing the sphere into the other two spheres. An example of this geometry for s = 2.1 is shown in Figure 7. The following quantities are measured: • Velocity U 1 of sphere 1, in direction of force.
• Velocity U 2 of spheres 2 and 3, in direction of force.
• Velocity U 3 of spheres 2 and 3, in direction perpendicular to force.

Wilson
Current method s  Table 1: Translational and angular velocities for triangle of spheres, computed using 64 × 64 grid, r/h = 3, κ = 5 and p = 20. Deviations from reference results by Wilson [35] are underlined.
As can be seen in Table 1, we are able to reproduce the results by Wilson to high accuracy for s ≥ 2.10 using a 64 × 64 grid. For s = 2.01 the results are inaccurate even though the quadrature is accurately resolved everywhere, and a finer grid is required if the results are to be improved. This is because the double layer density q gets sharply peaked as the particles get close to each other, with a peak that has spatial scale O(d 1/2 ) when the separation distance is d [6,32], and we can not hope to get good results if our grid is too coarse to resolve the density. The problem is therefore one of resolution, not quadrature. To solve the problem efficiently as s → 2, an adaptive refinement strategy for the surface grid would be required.

Periodic array of spheres
As a validation test for our periodic computations, we consider a periodic array of spheres of radius a arranged in a simple cubic lattice. This is modeled using a single sphere in a cubic box with sides L. We subject the sphere to an external force, and solve the mobility problem to get its velocity for various values of the concentration ρ = 4πa 3 /3L 3 . We then Figure 7: Flow velocity around triangle of spheres at configuration s = 2.1. Sphere 1 is the top sphere. A feature of boundary integral methods is that once the solution q is known, it is possible to zoom in and study the details of the flow field everywhere, as is done here in the right picture.
compute the dimensionless drag coefficient K = F/6πaU and compare it to the results available from Sangani & Acrivos [31]. Using a 30 × 30 grid and tuning all parameters (QBX+Ewald) for a tolerance of = 10 −7 , we get the results shown in Figure 8. Our results are identical to the reference results in all the digits presented in [31], all the way up to the limit ρ = ρ max , where the spheres touch (L = 2a). Our results for K also agree in all six decimals with the series expansion in χ = (ρ/ρ max ) 1/3 published in [31, eq. 60], up to and including χ = 0.5, after which the series expansion loses precision.

Application: Porous media flow
As a demonstration of a possible application for our method, we consider a periodic cell of dimensions 40 × 15 × 15 containing 131 randomly positioned, oblate particles with size (a, c) = (2, 1), for a volume particle concentration of about 24% system can be used as a model of a three-dimensional porous medium with non-spherical grains. To this end, we set a background flow u bg = (1, 0, 0) and solve the resistance problem for fixed particle positions. We use a 30 × 60 grid and tune our parameters for a tolerance = 10 −4 . The resulting system has 707,400 unknowns and is solved in about 2 hours on a regular workstation with a 3.4 GHz quad-core Intel Haswell CPU.
Once we have the solution q on the particles, it easy to evaluate the velocity field anywhere in the domain. We can then place a Lagrangian point randomly in the domain and let it be advected by the flow, using a high-order ODE solver (MATLAB's ode45).    [31]. The dimensionless parameter χ = (ρ/ρ max ) 1/3 measures how close the configuration is to the touching configuration (χ = 1).
Doing this for a few particles as they pass through the domain allows us to draw the streamlines of Figure 9a. Doing this for many (∼ 100) particles over a long time interval (t ∈ [0, 10 3 ]), it is possible to extract statistics about the flow. We shall here perform an analysis that closely follows that by de Anna et al. [13] for a two-dimensional medium with circular grains. Our intention here is not to draw any new conclusions about this particular case, only to show that our method makes the analysis possible and that the results correspond to those in [13]. The hydrodynamic particle dispersion can be measured through the particle dis- 2 , we can study the longitudinal displacement in σ 2 x and the transversal displacement in σ 2 yz := σ 2 y + σ 2 z . These are shown in Figure 9b, with time normalized by t A = d/v, the mean particle diameter over the mean velocity. Both σ x and σ yz can be seen to have an initial linear growth, until entering a transition regime beginning at t = t A , after which the growth approaches a rate proportional to t 1/2 (this is known as the Fickian regime). This is consistent with [13]. τ /tA = 4.00 Gaussian Figure 10: Probability density function of the normalized longitudinal Lagrangian velocity increments ∆ τ u/σ ∆u (τ ) when τ < t A and τ > t A . Clearly, a transition occurs around For a given time lag τ , the longitudinal Lagrangian velocity increment is defined as ∆ τ u = u(t + τ ) − u(t). The probability density function of the velocity increment, normalized by its standard deviation, is shown in Figure 10 for two different time lags. For large time lags (τ t A ) the distribution approaches the Gaussian distribution, while for small time lags (τ < t A ) the distribution is similar to that of the correlated continuous time random walk (CCTRW) model in de Anna [13].

Conclusions
We have in this paper presented a complete boundary integral framework for simulating periodic Stokes flow around spheroidal bodies. It is based on representing the flow using the Stokes double layer potential, which results in a well-conditioned system that converges rapidly when solving using GMRES. Singular and nearly singular integrals are computed using QBX, which we have adapted for the double layer potential. This allows us to evaluate the velocity to high accuracy everywhere in domain, such that we can have accurate quadrature also for nearly touching geometries. By exploiting the symmetries of the spheroidal bodies, we have been able to develop a precomputation scheme that minimizes storage requirements, while allowing us to rapidly compute the potential using QBX. By integrating the precomputed QBX with a fast and spectrally accurate Ewald summation method, we get a fast method which scales favorably as we go to systems with many particles.
The results presented in section 6 show that our method is able to accurately solve reference problems available in the literature. The example involving a porous media model also indicates that our method could be useful when studying such models, with both spherical and non-spherical grains. A natural application extension of our method would be to use it also for moving rigid particles, for example by studying sedimentation in particle suspensions, as was done in [2]. Then particles can come arbitrary close to each other, so close that lubrication forces start to dominate. This is not an issue for our quadrature, but may cause the double layer density to be sharply peaked, as we saw in the example in section 6.1.1. To deal with such cases, one could either (i) introduce an adaptive surface quadrature, that is able to resolve the density peak using a limited number of surface nodes, or (ii) couple the method to a lubrication model that removes the need for the double layer potential to resolve the strong interactions [21,24,25,32].
All density values are stored in the vector Q ∈ R N , The discrete approximation of the layer potential is then When solving a boundary integral equation, we typically need to compute the potential at all the points x i on the surface, which we can write more compactly as where A ∈ R N ×N is a (typically) full matrix and U ∈ R N is the layer potential at all the points on the surface, U i = u(x i ).

A.1.1 Rotational symmetry
Now, since the spheroid is rotationally symmetric and we use the trapezoidal rule with uniform spacing in the azimuthal direction, the grid will look identical from the point of view of all points which lie on the same latitude (θ). If we denote the rows of A by R i , i = 1, . . . , N , then this must mean that rows i and i + n θ contains the same information, and differ only by a permutation. This can be viewed as R i and R i+n θ being stencils that cover the spheroidal grid in the θ-ϕ-plane. Having their target point at the same θ-coordinate they then have identical coefficients, but they differ by a permutation due to the periodic wrap-around at φ = 2π. We introduce the permutation τ n (Q) 3 , such that We represent this permutation using a (sparse) N × N matrix P n , τ n (Q) = P n Q.
For 1 ≤ i ≤ n θ we then have that . . . , which we can summarize as So we only need to store the first n θ rows of A, as the remaining ones are given by successive permutations of that set. This reduces our memory use from N × N to n θ × N .

A.1.2 Mirror symmetry
For the rotational symmetry we could show that we only need to store coefficients related to the first n θ points in order to make computations for all n θ × n ϕ points. An additional reduction can be made by observing that the spheroidal grid also has a mirror symmetry around its equator. This means that we can use coefficients related to target points on the northern hemisphere to make computations for target points on the southern hemisphere. Hence, we only need to store the first n θ /2 rows of the matrix A in memory to be able to compute U from Q. Mathematically this is done for the current ordering of points by introducing a mirroring permutation matrix F n θ such that We can then recover rows n θ /2 + 1, . . . , n θ by using rows 1, . . . , n θ /2, and from there recover all the rows of the matrix using the rotational symmetry. For notational simplicity we here assume n θ to be even, extending to odd n θ is trivial.

A.2 Precomputing QBX
Now we have established the symmetries available when evaluating a layer potential on a spheroidal grid. We will now see how those can be used in conjunction with QBX, to create a fast method.

A.2.1 Laplace single layer potential
We begin by considering the Laplace single layer potential, having kernel K(r) = |r| −1 , as the analysis is simplified by all the involved quantities being scalar. Using the Laplace expansion (31) about an expansion center c, the potential can be written as where the expansion coefficients are computed as When this is done discretely on a spheroidal grid, there is a series of steps involved. The first step is to create a fine grid by increasing the number of grid points by a factor κ in each direction, so that we get a grid with points x i , i = 1, . . . , κ 2 N . The density is upsampled to the fine grid by a suitable interpolation scheme, which is represented as the matrix operation In the second step, we use the upsampled density to compute the local expansion coefficients up to order p at expansions points c i (located at a distance r in the normal direction from x i ), This R i has the symmetries discussed in section A.1. From rotation and mirroring when then have such that it is enough to directly compute R i , . . . , R n θ /2 , as the remaining R i can be recovered through these operations. The symmetries that apply to the rows R i also apply to the matrices M i , such that we can compute the coefficients at all expansion centers using M i for the first n θ expansion centers. Here, however, it is not enough to permute the right hand side vector; we also need to account for the phase of the factor e imϕ in the spherical harmonic. An expansion center c i+n θ will "see" all the grid points in the same angles ϕ as an expansion center c i does, plus a rotation of ∆ϕ. To account for this rotation, we introduce a phase matrix E ∈ C Np×Np with diagonal elements For 1 ≤ i ≤ n θ this lets us write . . . , and . . . , which can be summarized as From this we can recover the original rotational symmetry of R i , since E * E = EE * = I.
To also make use of the mirror symmetry for M i , we have to keep in mind that the mirroring operation makes an expansion center c n θ /2+i "see" the other grid points in a coordinate system which has been rotated 180 • about the x-axis, compared to what c i sees. This amounts to θ and ϕ in (95) having counterparts θ = −θ and ϕ = −ϕ in the mirrored system. Inserting this into the relevant terms in (95), we see that Changing sign in the exponential is equal to taking the complex conjugate of the entire expression. To handle the change of sign in the argument of the associated Legendre polynomial, we apply the parity relation available in Arfken & Weber [4, p. 776], Defining the diagonal matrix Θ ∈ R Np×Np such that where l = l(j) and m = m(j) as j = 1, . . . , N p , we finally arrive at the relation which allows us to compute all expansion coefficients for all expansion points using precomputed results for the first n θ /2 points.

A.2.2 Dipole potential
We will now extend the above results for the Laplace single layer potential to the dipole potential, defined in (32) for a vector density q as L [Γ, q] (x) := Γ q(y) · ∇ y 1 |x − y| dS y .
A discrete quadrature by expansion for this potential is formed from the expressions (33) and (34). We store the pointwise vector densities in linear form in Q ∈ R 3N , Just as for the single layer potential, the expansion coefficients are computed using an upsampled density Q, which here is computed as where ⊗ is the Kronecker product. The expansion coefficients are here computed as z lm (c) = κ 2 N n=1 w n q( x n ) · ∇ x 1 r l+1 Y m l (θ n , ϕ n ).
This can be represented in matrix form as where The difference from the single layer case is that the matrix S i ∈ C Np×3κ 2 N here includes the gradient of the spherical harmonic on the surface. Evaluation of the expansion is identical to the single layer case (103), with s(x i ) defined as in (100), and R i given by Due to the density now being vector-valued, there is difference between the single layer potential and the dipole potential in how we can make use of the grid symmetries.
To reuse M i of the first n θ points by rotational symmetry we must now not only take the permutation P n θ into account, but also the rotation of the frame of reference, since Q now represents a quantity which is pointwise vector-valued. Introducing the rotation matrix T z , the rotation of M i can be formulated as Similarly, for the mirror symmetry we must take the 180 • rotation about the x-axis into account. Representing this by the rotation matrix we arrive at the expression for the mirror symmetry. These results also carry over to the symmetry relations for R i ,

A.2.3 Gradient of expansion
To compute the gradient of a (single or double layer) potential that is locally represented by a local expansion at c, we differentiate (103) or (132) with respect to x, where g(x) ∈ C 3×Np has elements We can also write where g i (x) ∈ C Np , If we want to reuse g(x i ) from the first n θ points for the remainder of the points, we have to take into account that the output is vector-valued and has to be rotated into the correct frame of reference. So, if we have the entire process condensed into the form then ∇u(x i+αn θ ) = T α z (∆ϕ)R i P α n θ q, 1 ≤ α < n ϕ .

A.2.4 Stokes double layer potential
We now turn to our potential of interest, the Stokes double layer potential (8), Our goal is now to derive the matrix representations necessary for precomputing QBX for this potential on a spheroidal grid. We have from (38) that the we can represent the double layer potential as a combination of dipole potentials, which means that we can derive the expressions needed for the double layer potential using the results derived so far in sections A.2.1-A.2.3. Taking (148) apart, we get where u j (x) = L [Γ, q j n + n j q] (x) , j = 1, . . . , 3, u 4 (x) = L [Γ, y k q k n + y k n k q] (x) .
When computing these components on a spheroidal grid, we can improve the accuracy by taking into account that the geometrical quantities n and y are explicitly known, such that we only need to upsample the density q. Using (129), the expansion coefficients at a center c for the components u 1 -u 4 are then given by where the geometrical quantities on the fine grid are given by and P = diag (x 1 ·ñ 1 , . . . ,x κ 2 N ·ñ κ 2 N ) .
From the coefficient vectors z 1 -z 4 at c, the double layer potential at a point x can be computed as where the evaluation vectors s(x) and g i (x) are as defined in (100) and (143). For the points x i and centers c i related to the spheroidal grid, the operations (152), (153) and (158) can be represented using a set of transfer matrices M 1−4