{'fraction_non_alphanumeric': 0.04187313359367309, 'fraction_numerical': 0.01276539887769501, 'mean_word_length': 4.695010278452625, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient Implementation of Elastohydrodynamics via Integral Operators

The dynamics of geometrically non-linear flexible filaments play an important role in a host of biological processes, from flagella-driven cell transport to the polymeric structure of complex fluids. Such problems have historically been computationally expensive due to numerical stiffness associated with the inextensibility constraint, as well as the often non-trivial boundary conditions on the governing high-order PDEs. Formulating the problem for the evolving shape of a filament via an integral equation in the tangent angle has recently been found to greatly alleviate this numerical stiffness. The contribution of the present manuscript is to enable the simulation of non-local interactions of multiple filaments in a computationally efficient manner using the method of regularized stokeslets within this framework. The proposed method is benchmarked against a non-local bead and link model, and recent code utilizing a local drag velocity law. Systems of multiple filaments (1) in a background fluid flow, (2) under a constant body force, and (3) undergoing active self-motility are modeled efficiently. Buckling instabilities are analyzed by examining the evolving filament curvature, as well as by coarse-graining the body frame tangent angles using a Chebyshev approximation for various choices of the relevant non-dimensional parameters. From these experiments, insight is gained into how filament-filament interactions can promote buckling, and further reveal the complex fluid dynamics resulting from arrays of these interacting fibers. By examining active moment-driven filaments, we investigate the speed of worm- and sperm-like swimmers for different governing parameters. The MATLAB(R) implementation is made available as an open-source library, enabling flexible extension for alternate discretizations and different surrounding flows.


I. INTRODUCTION
Flexible filaments are ubiquitous in the natural world, and thus a clear understanding of their behavior is paramount in many biological problems. Models for simulating the dynamics of these filaments, while plentiful, have historically been mathematically complex and numerically expensive, with even simple computational experiments taking in the order of hours or even days to solve (see Moreau et al. for detailed benchmarking [1]).
Micro-scale filament problems have been previously tackled using other modeling approaches which can be broadly separated into those based upon (a) a nonlinear PDE in the filament position (such as the method by Schoeller et al. [2]), or (b) a discretization into simpler elements, such as beads with connecting springs [3] or interlocking gears [4]. For category (a), inextensibility is enforced using Lagrange multipliers of tension, which are often costly to compute. For the discrete approaches in category (b), other ways of enforcing this condition are used. For example, the bead model of Jayaraman et al. [3] prescribes large spring constants between each bead, contributing to the numerical stiffness of the system. Equivalently, the gears model of Delmotte et al. [4] imposes a non-holonomic constraint to ensure non-penetrability between adjacent beads, but as a result requires large numbers of points to represent a single filament.
A recent promising development via Moreau et al. [1], referred to as coarse graining, is based on reformulating the problem via an integral equation with the filament tangent angle as the dependent variable. The method, initially developed using a local hydrodynamic drag law, provides an efficient framework for simulating non-interacting filament dynamics. This approach builds upon the early studies of Brokaw [5,6] and Hines & Blum [7], and contrasts with Cartesian formulations [8,9].
The contribution of the present manuscript is to enable efficient and accurate simulation of multiple, non-locally interacting, passive and active filaments in ambient flows by incorporating recent developments in the regularized stokeslet method [10,11] with the integral formulation in terms of the tangent angle of Moreau et al. [1].
The potential applications for a fast and accurate filament modeling framework are numerous. There has long been interest in understanding the mechanics and regulation of sperm flagellar movement, in particular problems relating to: understanding the mechanical struc-ture and motor regulation [5,[12][13][14], investigating the response of the flagellar beat to its rheological environment [15][16][17], understanding the dynamics of sperm due to surrounding solid walls [18,19], and studying the effect of viscosity on sperm swimming [20]. For a detailed review surrounding the importance of the sperm flagellum see Gallagher et al. [21].
Furthermore, such a method could be used to investigate phenomena associated with epithelial cilia driven flows such as: cilia waveform modulation by length [5], the effects of flow induced by cilia on embryonic development [22], studying the physical limits of flow sensing [23], and investigating the mechanical structure of the axoneme in cilia [24]. Another avenue of active-filament research to which the proposed framework could be applied is magnetic swimmers [25]. These models have wider relevance in the field of synthetic biology, with particular application to microscopic bacteriophage-based fibre sensors [26][27][28] and flexible filament microbots [29]. The proposed framework could be used to further investigate the dynamics of bundles of filaments [30], and additionally has applications in the multi-scale studies of complex polymeric fluids, and flagellar movement through them [31,32].
We will extend the framework introduced by Moreau et al. [1], augmenting and reformatting their formulation with the method of regularized stokeslets of Cortez et al. [33,34]. These methods have proven to be accurate and effective in modeling the hydrodynamics in various multiple-fiber scenarios [35,36]. The method of regularized stokeslets enables the modeling of non-local hydrodynamics within and between filaments, and between filaments and surrounding structures. The method is implemented in a numerically efficient manner, retaining the computational economy and low hardware requirements inherited from the Moreau et al.

formulation.
The structure of this manuscript is as follows: in Sec. II, the Elastohydrodynamic Integral Formulation (EIF) for a single filament is proposed. In Sec. III, alterations to the EIF for various single-and multi-filament scenarios are presented. Verification and benchmarking of the method is given in Sec. IV. Simulation results of the problems formulated in Sec. III are then presented in Sec. V, followed by discussion of the results and of further possible applications in Sec. VI. The MATLAB R code for the methods described within this report are provided in the associated GitLab repository, available at gitlab.com/atticushm/eif.
The dynamics of elastic filaments in Stokes flow will be modeled by constructing integral operator formulations of the governing fluid and elastodynamic equations. Each filament is described by the time-evolving tangent angle θ(s, t) along its arclength. Taking X 0 (t) to be the leading point at time t, the filament geometry is then given by where s ∈ [0, 1] and θ ∈ [−π, π) are the dimensionless arclength (scaled by the filament length L) and tangent angle respectively. The velocity of a point on the filament is given by differentiating Eq. (1) with respect to dimensionless time, scaled by τ . The behavior of planar filaments in a Newtonian fluid is considered, which can be described by the three-dimensional dimensionless Stokes flow equations, where u(x, t) is fluid velocity, p(x, t) is pressure, and F (x, t) is a force exerted by the body onto the fluid. As shown by Cortez [34], an approximate solution to Eq. (2) is given by the regularized stokeslet integral, where f (s, t) is the force per unit length exerted by the filament on the fluid, nondimensionalized with the scaling µL/τ for a given fluid dynamic viscosity µ, and X(s, t) denotes the filament position as a function of arclength and time. The error arises from the regularization of the stokeslet, for a chosen regularization parameter 0 < ε 1, where r = 1 or 2 in the near-and far-field respectively. The combined process of solving Eq. (3) for the unknown force densities f (s, t) will result in errors of order ε everywhere [37]. Summation convention dictates that repeated indices are summed over and unrepeated indices range over {1, 2, 3}. The kernel of the integral in Eq. (3) is known as the regularized stokeslet, defined Applying the no-slip condition u(X(s, t), t) = ∂ t X(s, t) on the boundary of the filament yields the regularized stokeslet integral equation Inertialess dynamics requires that the filament is force and moment free at each instant, The integral formulation for the hydrodynamic problem thus comprises Eqs. (5), (6) and (7).
Next, the elastic behavior of the filaments is considered.
The integral operator for the elastodynamic behavior is formulated by considering the force and moment balance over an infinitesimal segment of a filament. We denote by N (s, t) and M (s, t) the contact force and contact moment respectively exerted by a distal segment of filament [s, 1] on a proximal segment [0, s). Constitutively linear elasticity and bending in the xy-plane implies that the dimensional contact moment relates to the curvature via where EI is the bending modulus of the filament. The fluid dynamic force density f (s, t) and N (s, t) are related by Under the assumption of free boundary conditions, the contact force and moment are zero at each end of a filament. Moment balance over the infinitesimal segment reveals that the contact force and moment are related as Integration of Eq. (10) over a distal segment [s, 1] yields, on application of the moment-free boundary condition and Eq. (9), Application of the the force-free condition N (1, t) = 0, the constitutive Eq.

A. Spatial discretization of the EIF
To solve the EIF, we spatially discretize filaments to obtain a set of integro-differential equations which can be numerically evaluated. Dividing the arclength domain into Q segments of equal length ∆s = 1/Q, the positions of the resulting segment endpoints are denoted as The angle connecting X [n] (t) to X [n+1] (t) approximating θ (n − 1)∆s, t is denotedθ [n] (t), for n = 1, . . . , Q. The positions of the endpoints are given in terms of the initial point X [1] (t) and discretized tangent angleθ [n] (t) as with the segment midpoints for each m = 1, . . . , Q. An illustration of this discretization is displayed in Fig. 1 For the fluid dynamics, rather than using a conventional (and potentially expensive), quadrature rule for evaluating the rapidly-varying kernel S ε jk , we employ the method of Smith [38]. By approximating the force density in Eq. (13) as piecewise constant along each segment, the kernel can be analytically integrated, reducing the level of quadrature needed to evaluate the slowly-varying force density (for higher-order force discretizations see Cortez [10]). Writing wheres [m] denotes the arclength at the midpoint of the m th segment, with a piecewise linear discretization for the filament, we obtain the spatially discrete equation, The integral in Eq. (24) can be calculated exactly by transforming to a local coordinate system in which one axis is aligned with the segment tangent (details are given in App. B, Eqs. B1-B3 of [38], although note that the simplified form for the near-field integral in Eq. B4 contains a typographical error in the δs/ε fraction, which is upside-down with the force and moment balance equations given in Eqs. (6) and (7) discretized via the midpoint rule as The semi-discrete form of the EIF is then ∂ tX [m] where and where m = 1, . . . , Q − 1 in Eq. (30) and m = 1, . . . , Q in Eq. (31) and (33).
By using Eqs. (32) and (33), the variablesX [1] (t), . . . ,X [Q+1] (t) and X [2] (t), . . . , can be eliminated from Eqs. (29), (30) and (31). The resulting system is then linear in Hence, given the discrete configuratioñ X [1] (t),θ [1] (t), . . . ,θ [Q] (t), the rate of change of position and angle can be found by solving a dense (3Q + 2) × (3Q + 2) system of linear equations. Thus, the semi-discrete system can be expressed concisely as an autonomous non-linear initial value problem, where Y (t) := [X [1] (t),θ [1] (t), . . . ,θ [Q] (t)] is a Q + 2 column vector andẎ (t) = ∂ t Y (t). By augmenting the problem with the unknown force densities, this can be written as the matrix system,  28)). The remaining zero entries correspond to the equivalence between the hydrodynamic and kinematic velocities (Eqs. (31) and (32)). The matrix system given in Eq. (35) is solved forẎ andf at each time step using the MATLAB R backslash command, and the resulting rates vectorẎ is integrated using the built-in variable-step, variable-order ODE solver ode15s [39]. To demonstrate the ease of application of the EIF framework, all simulations are performed using the MATLAB R R2019a default settings for ode15s. In particular, the absolute and relative error tolerances are 10 −6 and 10 −3 respectively. While this IVP exhibits some stiffness, it is less stiff than the systems produced by other methods, as the integral formulation avoid the need of additional Lagrange multipliers to ensure filament inextensibility. At each time step, the filament is constructed according to Eq. (32) with ∆s = 1/Q, ensuring that filament arclength is preserved over the course of the simulation.

III. SINGLE AND MULTI-FILAMENT PROBLEMS
In the following section, we apply the EIF framework described in Sec. II to problems involving single or multiple filaments in the presence of body forces, surrounding flow, or undergoing self-powered propulsion.

A. A single passive filament in shear flow
We investigate the dynamics of a single passive filament in a linear shear flow, Additionally, from the dimensional version of Eq. (16), together with the scalings defined above, we obtain the non-dimensionalized elastohydrodynamic integral equation where the dimensionless viscous-elastic parameter quantifies the ratio of viscous to elastic forces on a shear timescale, where EI is the bending rigidity of the filament. The apparent flexibility of a filament is completely characterized using V, with large values describing floppy fibers, and small values stiff fibers. The similarity between Eqs. (16) and (40) results in a similarly discretized form as in Eq. (30).
The hydrodynamic shear flow equation (Eq. (39)) is semi-discretized following the steps in The system of equations for a single passive filament in shear flow is thus given by Eq. (35) but with the alterations and where the vector of unknowns [Ẏ ,f ] T remains unchanged from Eq. (35).

B. A single passive filament sedimenting under gravity
The EIF method also allows for implementation of a body force to the system. In this section, the simulation of passive filaments sedimenting under gravity is considered. Assuming uniform mass per unit length ρ, the force per unit length due to gravity acting upon the filament is −ρge 2 . Non-dimensionalizing with respect to the filament length L, time scaling τ = µL 4 /EI, and force density µL/τ , addition of the gravitational force density produces the force and moment balance equations where X c (t) is the center of gravity of the filament, f (s, t) is the force per unit length the filament exerts on the fluid, and G is the elasto-gravitational parameter which, along with Eqs. (45) and (46), has spatially discretized form As in Sec. II A, we form a matrix system encoding Eqs. (49)-(51) along with Eqs. (31) and which can be solved to find the filament velocities ∂ t X 0 (t), ∂ tθ [m] (t) and force densities

C. A single active filament
For the problems considered in Secs. III A and III B, we can consider active filaments by including a time-dependent moment density term added to the elastodynamic formulation given in Eq. (12). The moment balance in Eq. (10) is extended where m(s, t) is the moment per unit length that drives the actuation of the filament. Continuing the non-dimensionalization as in Sec. II, with time scaling τ = ω −1 , yields the nondimensionalized elastohydrodynamic equation for a single active filament where S is the dimensionless swimming number, defined with ω being the angular velocity of the swimming beat prescribed to the filament. Note that the swimming parameter is different from the commonly-used sperm number Sp, which has dependence on a chosen resistance coefficient. Following work by Gadêlha et al. [8] and Montenegro-Johnson et al. [18], the active moment m(s, t) is set as a traveling wave with amplitude m 0 , wave number k and angular frequency ω. The semi-discretized form of for m = 1, . . . , Q. In the following section, we consider the how the presented framework can be extended to model systems of multiple filaments. with fewer than half of the degrees of freedom required by the BLM (see Fig. 2b).

V. RESULTS
In the following section, results from the numerical experiments outlined in Sec. III are pre- We initialize our filament as a straight rod with a small perturbation following the method of Young [41], writing the initial angular configuration as for initial angle θ 0 and small perturbation parameter ∆θ 0 . As discussed by Tornberg & Shelley [9], prescribing a small perturbation to a straight filament can drastically change the dynamics from rotational Jeffery orbits to interesting buckling phenomena.
In Fig. 3, we demonstrate how changing the value of V affects the amount of buckling a filament experiences. The fiber is perturbed with θ 0 = 0.9π and ∆θ 0 = 0.1. As the filament rotates, buckling modes form, which are directly linked to the size of V and the choice of perturbation (Eq. (58)). Large values produce high-order buckling modes, which can be seen in the fourth row of Fig. 3. Conversely, comparatively small values produce no buckling and the filament rotates through a standard Jeffery orbit. For V = 5 × 10 3 , a first-order buckling mode begins to form, visible in the first row of Fig. 3. The amplitude of the buckling increases as V increases until higher-order modes are induced. Tornberg & Shelley [9] examined buckling governed by an effective viscosity parameterμ = 8πV, producing filament shapes akin to the second row of Fig. 3.
The buckling problem of a flexible filament in shear flow has multiple solutions [40,42]. As noted by Tornberg & Shelley [9], the choice of initial perturbation to the filament shape can preferentially lead to one of the solutions. For example, the particular initial condition considered in [9] produces a reflected filament shape when changing the sign of the perturbation.
Additionally, for the high-order buckling modes, increasing values of Q must be chosen in order to resolve the filament dynamics. However, given an initial condition that does not uniquely determine a specific solution branch (as for this problem with large values of V), the buckled solution becomes sensitive to the choice of discretization. Details of the convergence of the method for this problem is contained within Sec. SI.1 of the Supplemental Material.
At critical values of V, unstable buckling can occur. For example, when V = 3.4 × 10 4 , a higher-order buckling mode initially forms (row three of Fig. 3) but collapses into a lowerorder mode. In Fig. 4, the filament evolution highlighted in Fig. 3 is displayed as a curvature plot in arclength and time. The unstable higher-order mode can be clearly seen collapsing into a lower-order configuration, indicated by the two dotted lines in Fig. 4a. Increasing V further ensures that a higher-order mode is stable throughout the compression phase (Fig. 4b).
The perturbation to the filament shape induced by buckling can be investigated by considering the evolving body-frame tangent anglesθ is the mean filament angle. By fitting Chebyshev polynomials toθ By sampling θ from the very-low amplitude parabola y = 1 × 10 −7 x 2 , the filament geometry is initialized with unit arclength, and pre-solved with a coarse discretization (Q = 10) until the shape is sufficiently curved so that a higher-resolution representation can be employed Streamlines indicate the direction of the surrounding fluid flow. In each case, the filament is modeled using Q = 40 segments, with initial shape parameters θ 0 = 0.9π and ∆θ 0 = 0.1 (Eq. (58)). For a given initial condition, the size of V completely determines the level of buckling the filament exhibits as it is compressed during the rotation. For large values, high-order buckling modes appear, as in the third and fourth rows of the above diagram.  For different choices of G, various sedimentary buckling modes can be observed. In Fig. 6, a stable "U" filament configuration emerges over time. For large values of G > 3000, a metastable "W" configuration forms, which transitions into a stable "U" horseshoe shape (see Fig. 7, in which G = 3500). Such behavior has previously been observed by Delmotte et al. [4] and Cosentino Lagomarsino et al. [43], who witnessed buckling at the same threshold value for their identical elasto-gravitational parameter. This transition shifts the filament's center of gravity to the left, creating an asymmetry which is then partially resolved as the horseshoe equilibrium configuration is reached.

C. Results: a single active filament
In the following section, we consider swimming in a stationary fluid caused by two types of traveling-wave moment densities, (1) sperm-like sinusoidal motion m 1 (s, t), and (2) a wormlike beating pattern m 2 (s, t), given respectively by in which T = 2π is the period of the driving wave and X Swimming speed for different choices of parameter pairs (S, m 0 ) are presented in Fig. 8. For critical values, filaments self-intersect, in which case the EIF is inapplicable. The shape of such filaments are shown in Fig. 9. For a sperm-like swimmer (left panel of Fig. 8), swimming speed increases as m 0 and S are increased in tandem for a fixed wave-number k = 4π. This is in contrast to worm-like swimmers (right panel of Fig. 8), in which there is a clear optimal choice of m 0 for a given S to induce fastest swimming. Shape profiles for the fastest swimmers of each swimming type are also displayed in Fig. 9. The shape difference between two filaments is quantified using the Procrustes measure, calculated using the procrustes MATLAB R command. Higher values indicate a larger degree of dissimilarity between the two filaments [44].

D. Results: multiple passive filaments in shear flow
Motivated by the work of Young [41], two filaments of equal length are placed into a linear shear flow so that their initial midpoints X(s = 1/2, t = 0) intersect the line y = 0, and are separated by a distance ∆X 0 . The initial shape profiles for each filament are the same as in Sec. V A. For two filaments characterized by the same viscous-elastic parameter, non-local hydrodynamic interactions lead to geometric asymmetry, as seen in Fig. 10. The value of V determines the level of dissimilarity between filaments with identical initial shape profiles, with higher values leading to larger deviations in shape throughout the rotation.
The initial separation distance also has a significant influence on the dynamics of each filament. The shapes of initially close filaments evolve in tandem, assuming similar geometries at a given moment in time (row one of Fig. 10). Increasing the separation distance results in a decoupling of the filament shape profiles (row two of Fig. 10). The EIF allows for each filament's characteristic parameter to be independently chosen, and in Fig. 11, we highlight the variation in shape this can induce. In each of the panels in Figs The geometric similarity between filaments can be assessed via their Procrustes score [44], calculated using the procrustes MATLAB R command. We compare pairs of filaments initially separated by ∆X 0 = 0.5 and 1 for a range of V, with results displayed in Fig. 12. As steady-state profiles being reached sooner than in the single-filament cases. Furthermore, the onset of buckling occurs at reduced values of G when multiple filaments are interacting.
As was the case for multiple filaments in shear flow, the inter-filament spacing ∆X 0 also has a significant effect on the resulting group dynamics, and can lead to symmetry breaking in the arrangement of filaments. In Fig. 14, two initial filament placement configurations are tested.
For each setup, we initialize four filaments arranged in a grid, with ∆X 0 = 1.5 and vertical separations ∆Y 0 = 0.5 (row one of Fig. 14) and ∆Y 0 = 1.5 (row two of Fig. 14). Smaller vertical spacing leads to the filaments "nestling" in a horizontally-mirrored configuration; For all choices of G, vertical symmetry-breaking occurs between the top-most and middle rows of filaments, with the second and third rows of filaments on the left and right of the frame nestling into those below them. For G = 3500, more prominent buckling is apparent, with filaments in the central column approaching a "W" shape, whereas those on the flanks assume a "U" configuration. Competing interactions between filaments in the central column and their surrounding neighbors causes them to remain in a metastable "W" shape longer  16. Two identical filaments with S = 8 and m 0 = 0.05, propelled by a sperm-like active moment, swimming alongside one-another. In (a), filaments in close proximity of each other produce a higher average velocity (VAL). In (b), the shapes of the evolving filaments are presented for an initial vertical separation ∆Y 0 = 0.1 for t ∈ [6π, 24π]. As in Fig. 9, the magenta paths in (b) indicate the path the leading point of the filament takes over the course of the simulation.
than would be expected in the single-filament case.

F. Results: multiple active filaments
Multiple self-propulsive active filaments can be simulated using the EIF method. To illustrate this, we consider two filaments swimming alongside each other at different separation distances ∆Y 0 , with shapes initialized as in Sec. III C. For fixed S = 8 and m 0 = 0.05, the size of the separation determines the resulting average swimming velocity, as shown in Fig. 16.
At low levels of separation, hydrodynamic interactions between the synchronous beats result in a higher average VAL (Eq. (62)), which decays as the filaments are placed further apart (Fig. 16a).

VI. DISCUSSION
In this manuscript, the elastohydrodynamic integral equation (EIF) framework is formulated, and applied to problems involving single and multiple filaments in shear flows, sedimenting under gravity, and swimming due to a prescribed active moment. From the simulations presented in Sec. V, it is apparent how inter-filament non-local hydrodynamic interactions have a role in governing filament shapes and buckling behavior. By examining active moment driven swimmers, an optimum pairing of moment-amplitude and characteristic swimming number is revealed.
One of the key benefits to the integral formulation framework presented is that the need for computing Lagrange multipliers of tension is removed; inextensibility of the modeled filament is already guaranteed by the use of the method of lines discretization. This results in a method that is more computationally efficient than many similar contemporary models.
The EIF method proposed by Moreau and co-authors [1] highlighted the reduction in computational runtime attainable when using the integral formulation applied to elastodynamics, due in part to the reduction in the required degrees of freedom the integral formulation affords. Alleviating this numerical stiffness facilitates the study of non-linear problems, such as those involving filament buckling investigated in Sec. V. Additional computational costs are incurred by introducing non-local hydrodynamics to the EIF framework, in line with the increase in the dimension and density of the matrix system solved at each time-step. The BLM and the EIF with regularized stokeslets complete a relaxing rod simulation in about equivalent walltime (as expected, since both approaches compute three-dimensional non-local hydrodynamics), shown in Fig. 2c. However, to ensure that the point-wise error is small, fine discretizations need to be used with the bead model (Q > 60), whereas the EIF method performs equally well with a much coarser segmentation of the filament (Q > 20). This allows for adequately accurate results to be obtained for reduced computational cost when using the EIF with regularized stokeslets over the bead and link model.
The proposed method has the potential to quickly and accurately simulate arrays of filaments in various flows and surroundings. The EIF method presented here will enable the solution of more challenging problems such as planar beating cilia sheets, or multiple sperm swimming in a narrow channel. The methods' modular framework enables such problems to be setup and executed in a simple and obvious manner. Additionally, the EIF method fits into the family of regularized stokeslet methods [34], which are increasingly widely used for problems in biological fluid mechanics. Crucially, the computational efficiency exhibited in the single filament experiments is retained. The method could be extended by incorporating the treecode formulation of Wang et al. [45], or equivalent methods, to coarse-grain the farfield flow, simplifying the computational cost and enabling simulation of increasingly large numbers of filaments. The inclusion of a repulsive force term, such as the Lennard-Jones potential employed by Jarayaman et al. [3], could additionally enable simulation of filaments in close proximity and help avoid filament self-intersection.
The results of Sec. V A suggest that the shape of a buckling filament can be described to a high degree of accuracy by a relatively low order Chebyshev polynomial (Fig. 5). These results suggest that a modified discretization based on orthogonal polynomials, perhaps in concert with suitable quadrature techniques [46] might provide further improvements in efficiency and scalability. In summary we hope that the integral operator formulation of elastohydrodynamics will be valuable in biological fluid dynamics and beyond.