Theoretical framework for two-microswimmer hydrodynamic interactions

Hydrodynamic interactions are crucial for determining the cooperative behavior of microswimmers at low Reynolds numbers. Here we provide a comprehensive analysis of the scaling laws and the strength of the interactions in the case of a pair of three-sphere swimmers. Both stroke-based and force-based elastic microswimmers are analyzed using an analytic perturbative approach, focusing on passive and active interactions. The former are governed by the cycle-averaged flow field of a single swimmer, which is dipolar at long range. However, at intermediate distances, with a cross-over at the order of 102 swimmer lengths, the quadrupolar field dominates which, notably, yields an increase of the swimming velocity compared to individual swimmers, even when the swimmers are one behind another. Furthermore, we find that active rotations resulting from the interplay of the time-resolved swimming stroke and the ambient flow fields and, even more prominently, active translations are model-dependent. A mapping between the stroke-based and force-based swimmers is only possible for the low driving frequency regime where the characteristic time scale is smaller than the viscous one. Finally, we find that the long-term behavior of the swimmers, while sensitive to the initial relative positioning, does not depend on the pusher or puller nature of the swimmer. These results clearly indicate that the behavior of swarms will depend on the swimmer model, which was hitherto not well appreciated.


Introduction
Locomotion of microscopic organisms such as bacteria or sperms is governed by laws which are different from those governing the world tangible by humans. Due to the dominance of viscous drag over inertia and the ensuing time-independence of the Stokes equations, a successful swimming strategy has to break the time-reversal symmetry [1]. Several theoretical models [2][3][4][5][6][7][8][9][10][11][12][13], experimental realizations [14][15][16][17][18][19][20][21][22][23] and simulations [24][25][26][27][28] have been employed to scrutinize the details of locomotion under these laws, yet all these approaches can be subdivided into two main groups with respect to their driving. The first class, which has been investigated mostly, consists of swimmers with an inherently prescribed stroke, insensitive to the environment of the swimmer [4][5][6][7]. The second class comprises swimmers with prescribed forcing rather than prescribed stroke, where the swimmer's elasticity controls how strong the swimmer adapts its stroke to its environment [3,8,9,13]. However, at the level of a single swimmer typically a mapping between force-based (FB) and stroke-based (SB) swimmers exists, obtained by extracting the stroke in the FB model.
A self-propelled swimmer typically induces a dipolar flow field, allowing to classify swimmers as pushers or pullers. In addition to it, higher order flow fields are induced [29]. Interestingly, living microswimmers are capable of not only sensing these flows, but also reacting to them, which is associated with important physiological functions [30][31][32][33]. Trying to understand these functions as well as to explain the collective behavior of swimmers led to studies of hydrodynamic interactions between microswimmers. Here, the focus so far was on swimmers with prescribed stroke [34][35][36][37][38][39]. Interactions between two FB swimmers were, on the other hand, studied for effective microswimmer models [40][41][42], microswimmers without intrinsic elasticity [43] or hybrid SB-FB models [44]. A study for fully force-driven microswimmers, which necessitate an elastic degree of freedom, is yet pending. Consequently, it is still not understood whether microswimmer interaction is preserved under the mapping between FB and SB models.
Earlier work on SB swimmers established the distinction between the so-called passive and active interactions [34,37,39]. Accordingly, passive interactions refer to the contribution to the translational or angular velocity which would similarly be experienced by a passive, undriven object of the swimmer's average shape. Thus, passive interactions are independent of the swimmer's own activity and are solely due to the flow fields produced by other swimmers. In addition to this, a second component to the interactions has been reported, which explicitly depends on the swimmer's own activity [34,37,39] and its interplay with the ambient flow fields. This motivates the terminology of active interactions for the second contribution, although clearly both active and passive interactions root in the active nature of microswimmers. The total translational swimming velocity thus can be written as v tot. = v 0 + Δ v pass. + Δ v act. , ( 1 ) with v 0 the swimmer's self-propulsion speed in an otherwise empty fluid, and Δ v pass. and Δ v act. passive and active interactions. The angular velocity Ω decomposes similarly. The same concept is also applicable to FB microswimmers, however a meaningful definition of passive interactions in this framework does not refer to the swimmer's average shape, but to how the elastic undriven swimmer itself behaves in the flow field produced by other swimmers, with no additional constraint on its shape. The linear superposition principle of low Reynolds number hydrodynamics predicts that rigid swimmers which propel as a result of a prescribed surface slip velocity, e.g. so-called squirmers, react to ambient flow fields similarly as a passive swimmer does [36]. Squirmers, which are typically used to model ciliated microswimmers as Paramecium, or algae colonies such as Volvox, thus interact by definition only passively. However, in opposition to the expectation that passive interactions are in general a sufficient description, it has been reported for the linear three-sphere swimmer that the active component is actually dominant for swimmer separations smaller than a threshold value which depends on the details of the stroke [34]. Such active interaction effects therefore might be of significant relevance also for biological microorganisms that propel by shape changes, as e.g. Chlamydomonas or euglenoids, and therefore also might shed new light on their collective behavior. However, several competing models for the active interactions are currently discussed, providing different predictions for the scaling and the sign of the interaction effects [34,37,39,44]. Given that most models rely on the same basic assumptions, the origins of these differences have not been clearly established thus far.
We here address these issues using a pair of linear three-sphere swimmers. This design has proven to be particularly useful already at the level of a single entity. Its simple design allows for analytic tractability, yet the main physical principles of self-propulsion are correctly captured in both the SB and FB models. Therefore, the results obtained for these swimmers could be applied in very different contexts. Hoping that similar universality of trends can be also established for hydrodynamic interactions, we now study microswimmer interactions in both the FB and the SB models using a recently developed perturbative approach [9] (section 2). Details on the perturbative calculations for the FB and SB models as well as for the proportionality of the leading order swimming velocity and the time-averaged flow fields are given in the appendices A-D.
The above calculations allow us to resolve the similarities and differences between these two different families of swimmers (section 3). Notably, after calculating the average flow fields of our swimmers, we find equivalent passive interactions (Δ v pass. ) in the two models, within the leading orders in stroke amplitude and inverse swimmer distance considered here. To the contrary, the active component (Δ v act. ) depends on the swimming strokes of both devices and differs already at the leading quadrupolar order. Interestingly, the interactions between FB swimmers become equivalent to the interactions between SB swimmers if the driving frequency in the FB model is small as compared to the inverse viscous time of the system. Hence, the interactions of SB swimmers do actually present a special case of the interactions of two FB swimmers. These findings are systematically discussed in the context of the existing literature [34,37,39,44], clarifying the observed differences in the scaling laws of the interaction strength. This analysis is complemented by an in-depth discussion of the effect of the hydrodynamic coupling on the propulsion velocity. Unlike previously observed, we show that two collinear or side-by-side swimmers with parallel swimming direction typically benefit mutually from each other due to the time-averaged flow fields they produce. In simple words, a swimmer propagates faster in a pair compared to when it is alone, even if it is leading. This result, which is a consequence of vanishing Reynolds numbers, is independent of the details of the swimming stroke due to a close relation of the leading order swimming velocity and the time-averaged flow fields produced by each swimmer. The most important findings are summarized in section 4, which concludes the paper.

Force-based model of the three-bead swimmer
We consider two linear three-bead swimmers identified by subscripts s, p ∈ {I, II}. We restrict our analysis to a two-dimensional configuration space (d = 2) for the sake of simplicity, in our case the x-y-plane, as depicted in figure 1. Each swimmer consists of three identical spherical beads of radius a, denoted by subscripts i, j, m ∈ {1, 2, 3}. The beads are connected by identical linear harmonic springs of spring constant k and equilibrium length L, such that the spring force between two connected beads si and sj, acting on the latter, is given by with x = R sj − R si the vector connecting both beads. We employ a convenient double index notation to identify the beads, where the first index corresponds to the swimmer and the second index to the bead number within this swimmer. In this work we only consider swimmers with a rigid joint connecting both swimmer arms, i.e. the two swimmer arms cannot bend with respect to each other. This allows to systematically compare FB to linear SB swimmers [5,34,37,39], to which this constraint is inherent.
As shown in figure 1, swimmer I is oriented along n I = (sin θ I , − cos θ I ), swimmer II along n II = (sin θ II , − cos θ II ) and the middle beads of both swimmers are separated by r = r(sin(ϕ + θ I ), − cos(ϕ + θ I )). Here, ϕ denotes the angle enclosed between n I and r. The swimmers are immersed in a viscous fluid of viscosity η, where the Reynolds number of a single bead is assumed to be zero. The interaction between two spherical beads at sufficiently large distance is then described by the Oseen tensorT with 1 the unit matrix and ⊗ the tensor product. Each swimmer s is driven by sinusoidal forces of frequency ω acting on the two outer beads (bead s1 and s3) along the respective adjacent arm with relative phase shift α s . The driving forces on swimmer II precede those on swimmer I by a phase shift of γ: with E si (t) the force on the bead si, A si the respective force amplitude and t the time. We introduce a second set of independent variables for the driving forces which is required in the perturbative calculation, A the overall strength of the driving forces and B I , B II , C dimensionless relative parameters. The force on each middle bead s2 is determined by requesting the sum of all three forces on the beads of swimmer s to vanish [8], E s2 (t) := − E s1 (t) − E s3 (t), as a prerequisite for self-propulsion. We then devise an equation of motion for all of the n = 6 beads, with R si (t) the position of the respective bead and μ = (6πηa) −1 the Stokes mobility. NN(j) denotes all beads connected to bead j via a spring within the geometry of the linear swimmer, and in the remaining summation at the beginning of the second line we sum over all beads pj different from the bead si. To obtain the swimmer behavior analytically, we employ a recently developed perturbative approach to bead-spring microswimmers [9] (details in appendix B). Also, we solve (5) numerically using the NDSolve function from Mathematica [45], with a superimposed angle spring potential ensuring the rigidity of the joint at the middle bead of each swimmer (details in appendix A).
A system of bead-spring swimmers as described above can be rescaled using the bead radius a as the length unit and the viscous time t V := 6πηa/k as the time unit, in order to identify the effective parameters [9]. We define the dimensionless parameters q := a/r and ν := a/L, encoding together with ϕ, θ I and θ II the geometry of the system, and the rescaled driving frequency Γ := ωt V , comparing the two time scales set by the driving forces and the viscous time.
In the subsequent analysis, we will illustrate the results for the swimmer interactions using a swimmer of puller-and a swimmer of pusher-type with distinct parameters. For the puller, we choose and for the pusher-type swimmer Those rescaled parameters are derived from a set of parameters which optimizes the computational costs in lattice-Boltzmann simulations, while ensuring sufficient accuracy of the results. Similar parameters have been used in simulations of a single three-bead swimmer [46], and we adhere to them to simplify comparison to lattice-Boltzmann simulations of many swimmers in the future. The values for the driving force amplitudes, A s1 and A s3 , have been chosen in their order of magnitude such that they correspond to moderate arm oscillations of approximately one bead radius, maximizing the swimming speed, and are small enough to avoid the beads coming closer than 3-4 bead radii, where the Oseen approximation becomes invalid.

Stroke-based model of the three-bead swimmer
In the SB model, we consider two swimmers with identical geometry as in the FB model. In contrast to the FB swimmer, the length of each swimmer arm is directly prescribed as a function of time with L the average arm length, ξ sb the corresponding arm oscillation amplitude (indices for the swimmer arms b, c ∈ {1, 2}) and β I , β II and δ the phase shifts within each swimmer and between both, respectively.
Similarly to the FB model, we introduce a second set of parameters for the arm oscillation amplitudes required in the perturbative analysis, ξ := ξ I1 , D I := ξ I2 /ξ I1 , F := ξ II1 /ξ I1 , D II := ξ II2 /ξ II1 , with ξ the overall amplitude of the arm oscillations. With both swimmers constrained to linear shape and to the prescribed arm lengths, the system contains six undetermined degrees of freedom, which we choose to be the positions of both middle beads, R I2 and R II2 , as well as the orientations of both swimmers, θ I and θ II . The positions of all other beads are then given by the prescribed arm lengths and the constraint on the linear swimmer shape.
The relation of the bead velocities to the forces on the beads is given similarly as in the FB model by where F si (t) denotes the forces which are acting on the beads to enforce the prescribed strokes. The remaining six degrees of freedom are then determined by the condition that the total vectorial force on each swimmer vanishes, as well as the total torque on each swimmer, which is a scalar quantity in the 2D framework employed here: Using the compact notation R = ( R I1 , R I2 , R I3 , R II1 , R II2 , R II3 ) and F = ( F I1 , . . . , F II3 ), we can re-express (9) as d dt Here, μ(R) denotes the (n · d) × (n · d)-dimensional mobility matrix [9,47] and we omit the time-dependence of R and F in our notation for the sake of brevity. Vectors in the (n · d)-dimensional configuration space of all bead positions are denoted by bold symbols and higher order tensors on this space by underlined symbols. In order to perturbatively calculate the swimmer behavior, we expand the undetermined variables R I2 , R II2 , θ I and θ II as power series in ξ and q, and from this calculate an expansion of dR/dt and μ(R). By also expanding (12), we are able to express the components of F associated to different powers of ξ and q in terms of the middle bead positions and swimmer orientations. The force-free and torque-free conditions, (10) and (11), close the equations obtained and by solving them we obtain the full swimmer behavior as an expansion in ξ and q. The perturbative calculation is explained in more detail in appendix C.2.

Mapping from force-based to stroke-based model
For a single FB swimmer, the swimming velocity can be either calculated directly in the FB framework [9], or one can equivalently extract the time-dependent arm lengths and insert them in the respective expression for the swimming velocity in the SB model [48]. Since the single swimmer behavior is uniquely defined by the time-dependent arm lengths together with the force-free condition [48], both ways must yield the same result. This means that it exists a mapping from the parameter space of the FB system to the parameter space of the SB system which preserves the overall swimmer dynamics. For the sake of simplicity, we restrict the mapping to the first order in A and ξ in the subsequent calculations, where the average arm lengths in the FB model are still equal to the spring lengths in the mechanical equilibrium, L [9]. Hence, we assume the geometry of the FB swimmer in its mechanical equilibrium for the time-averaged shape of the SB swimmer, i.e. we assume average arm lengths L, and only require to map the swimming stroke, i.e. the arm oscillation amplitudes. This reproduces consistent second order swimming velocities and flow fields in both models and is hence sufficient for the subsequent comparison.
The inverse map from the SB to the FB system is not unique, since the FB model comprises one additional parameter, namely the spring constant k. It gives rise to the viscous time scale t V and the rescaled driving frequency Γ, with respect to which a FB swimmer attains the maximum driving speed at Γ ≈ 1 [27], and is unable to self-propel for Γ 1 and Γ 1 [49]. In contrast, ω has no impact on the propulsion of a SB swimmer measured over one stroke [48] and drops out in a suitable rescaling. This is because with a prescribed stroke, calculating the swimmer behavior reduces to a purely geometric problem due to the time-independence of the Stokes equations [48].
The natural question arising then is: is also the interaction of two or more swimmers equivalent in the FB and SB models under the above mapping, or does the interaction in the FB model change with Γ, when the driving forces on each swimmer are such that the stroke amplitude and phase of each swimmer alone would be constant?

Flow field of a single swimmer
A single linear three-sphere swimmer with equal bead radii typically produces a flow field which is dominantly dipolar in the very far field and quadrupolar at intermediate distances ( figure 2) [34]. The transition between both fields depends on the swimmer parameters, and appears in our case at the order of thousands of bead radii or equivalently hundreds of swimmer lengths. The time-averaged dipolar flow field scales to leading order as A 4 or ξ 4 , depending on the model considered, and is given in terms of the swimming stroke as with r := | r| and n the unit vector along the swimmer axis (details on the calculation in appendix C.1). We have omitted the swimmer index in ξ 1 and ξ 2 since this result holds for a single swimmer. In contrast, the average quadrupolar flow field scales to leading order as A 2 or ξ 2 and reads in terms of the stroke Using the mapping from the FB to the SB model at sufficient order in ξ and A, one obtains similar expressions for the FB swimmer, which we do not include as they would become too large.
To understand the predominance of the quadrupolar regime in the average flow field up to distances of hundreds of swimmer lengths, we first consider the special case of a swimmer with ξ 1 = ξ 2 , i.e. equal arm oscillation amplitudes. In this case, the swimmer becomes invariant under a combined time-reversal (t → −t) and parity ( r → − r) transformation (TP transformation) and the dipolar regime in the time-averaged flow field is lost [34]. Indeed, flow fields which decay with an even order in the distance r from the swimmer (dipolar, octopolar, etc) are not consistent with TP invariance. Thus, the time-averaged far field of a TP invariant and self-propelled swimmer is quadrupolar, decaying as r −3 , and the dipolar regime in figure 2 would disappear. However, even a TP-invariant three-sphere swimmer still produces a non-zero dipolar flow field that oscillates with the frequency of the swimming stroke.
For ξ 1 = ξ 2 , which is the case for both the puller-and pusher-type swimmer defined in section 2, the TP invariance is broken and a non-zero time-averaged dipolar flow field arises. Notably, the leading order, ξ 2 contribution to all even-order time-averaged flow fields still vanishes, and the first non-zero contribution arises at order ξ 4 . The reasoning to show this extends an argument made originally by Golestanian and Ajdari [48]. Similarly to the time-averaged swimming velocity, the amplitude of each component (dipolar, quadrupolar, etc) of the time-averaged ξ 2 flow field must be given by a geometric prefactor times the area in the configuration space enclosed by the swimmer's stroke (appendix D). First, this leads to the important conclusion that, at order ξ 2 , average flow field and swimming velocity are proportional and linked by a factor depending on the average swimmer geometry only. Accordingly, knowing the swimmer's geometry, one can, at leading order, directly infer the average flow field from the swimming velocity. Second, we conclude that the proportionality coefficient needs to vanish for all even-order components (dipolar, octopolar, etc) of the average flow field, since it vanishes for the special case ξ 1 = ξ 2 and is generic for linear three-bead swimmers with equal radii and equal average arm lengths.
In the case of sinusoidal driving forces with a single frequency, also the ξ 3 time-averaged flow field vanishes since it is associated to odd products of the sinusoidal stroke, which average to zero over a stroke cycle. This explains why passive dipolar interactions therefore arise firstly at fourth order in the stroke amplitude (∼ξ 4 /r 2 ) (13). In [34], a stroke comprising stepwise arm contractions at constant velocity [5] has been employed, associated to infinitely many harmonics of the base frequency ω in a Fourier decomposition. Therefore, a non-zero time-averaged dipolar flow field has already been observed at third order in ξ, making the quadrupolar regime less dominant than in the case of sinusoidal driving with a single pure frequency.

Passive interaction
While a single swimmer alone would propagate with velocity v 0 along its axis, in the presence of another swimmer it experiences both additional translation and rotation due to the hydrodynamic interactions. We split the swimmer's translational velocity into a component v 0 + Δv par parallel to its axis n s and into a component Δv ort along n ⊥ s (see figure 1). Additionally, the swimmer rotates with an angular velocity ΔΩ which is defined as positive when the swimmer rotates counterclockwise (CCW). We restrict our subsequent analysis to the velocities time-averaged over one swimming cycle.
We find that the interactions experienced by swimmer I, without restriction of generality, compose of a component which is independent of the swimming activity of swimmer I and the phase shift between both swimmers, and a second component which is not and does depend on the phase shift, γ or δ in the FB and SB models respectively, via a sine or cosine function. The first component is identified as passive interaction, according to the definition given in the introduction, while the second component clearly has to be active. We find that passive translational interaction is to leading order given by with R I2 the current time-averaged position of the second bead and u the time-averaged flow field produced by the second swimmer. The passive rotation experienced by a linear swimmer is to leading order given by with n s and n ⊥ s the current time-averaged orientation of the swimmer and the normal to it. This expression for the linear swimmer can be understood analogously to the vorticity of the flow field for a swimmer of spherical shape: while a spherical swimmer senses the gradient of the flow field in all spatial directions, the linear swimmer senses it only along its axis. Note that for the passive interactions experienced by swimmer I the single-swimmer flow field of swimmer II is relevant, and not the average flow field produced by swimmer II in the presence of swimmer I. Since the difference between the two situations arises from the swimming activity of swimmer I, this effect gives rise to active interactions.
Mediated by the passive interactions, the transition from predominantly quadrupolar to dipolar with increasing distance is also found in the interaction of two collinear swimmers, i.e. swimmers with a common axis (θ I = θ II = 0, ϕ = 0) (figure 3). For side-by-side swimmers with parallel swimming direction, defined by θ I = θ II = 0 and ϕ = π/2, passive dipolar interaction effects are found only in v ort and passive quadrupolar interactions only in v par (figure 4). This is due to the shapes of the respective flow fields (figure 2). In the general case of two interacting swimmers, also passive ξ 4 or A 4 quadrupolar rotation would be observed, resulting from the gradient of the dipolar time-averaged flow field. However, in the two most simple swimmer configurations chosen here for illustration, this term is absent due to the symmetries of the dipolar flow field.
The proportionality of the ξ 2 or A 2 time-averaged flow field and swimming velocity has important implications for the swimmer interactions at this order. Namely, for swimmers of given geometry, the passive interactions between the swimmers depend only on their swimming velocities and their relative positioning, but not on the details of the strokes. The factor of proportionality for the quadrupolar flow field of a single swimmer, measured in front of the swimmer along the swimmer axis n s , relative to the swimming velocity is given by This ratio is positive for values L 6a where the Oseen approximation can be expected to be sufficiently good, assuming that the bead's oscillation amplitude is not substantially larger than one bead radius a, i.e. A si /ka 1. Therefore, two linear swimmers in the collinear or side-by-side configuration with the same swimming direction always mutually benefit from passive interactions in the quadrupolar regime (figures 3 and 4). While in some of the previous works the same direct correspondence of passive interactions and time-averaged flow field produced by surrounding swimmers has been found [34,39], other authors have reported passive interactions different from the time-averaged flow field. In particular, for a similar swimmer model passive quadrupolar interactions scaling to leading order as L 1 for L → ∞ have been reported [37,44], while the time-averaged quadrupolar flow field produced by a single three-sphere swimmer scales to leading order as L 0 in the same limit [34]. We have been able to reproduce the linear scaling in L of the passive interaction in a slightly altered SB perturbative calculation, where the distance between the two swimmers is assumed to be a constantly r throughout one swimming stroke, instead of treating it as an unknown variable to be determined. In this work however, we impose no constraints on the distance between both swimmers, only assuming that the initial distance at t = 0 is r. Consequently, we obtain the exact correspondence of passive interaction and the average flow field.

Active translation
Beyond the approximation of the swimmer as a passive object drifting in the local time-averaged flow field around it, the swimmer deforms actively in order to self-propel, introducing additional active interaction effects which depend on the swimming stroke of both interacting swimmers as well as their relative phase. In both the FB and the SB model, active translational interactions arise firstly at the quadrupolar order (∼ 1/r 3 ) and order A 2 or ξ 2 , depending on the choice of model. In both models, this term has only a component parallel to the swimmer axis for arbitrary swimmer positioning. The analytical result for the change in swimming velocity of swimmer I in interaction with a second swimmer II side-by-side with it reads in the SB model including both active and passive terms (details of the SB perturbative calculation are presented in appendix C.2). For the FB model, we restrict the analytical result to two pullers in the same configuration, as otherwise the expression would become too large, For the details of the FB calculation we refer to appendix B. Both results, (18) and (19), scale to leading order constant in L in an expansion around L = ∞, similarly to the passive effect (14). For moderate stroke amplitudes of both swimmers (A ≈ ka, ξ ≈ a), we therefore observe that passive and active quadrupolar interaction are typically of the same order of magnitude. In the collinear and side-by-side configuration, passive interactions always enhance the self-propulsion of both swimmers, whereas the sign of the active interactions depends on the phase shift between both swimmers via a sine or cosine function, as shown in (18) and (19). Thus, for a large part of the parameter space, two swimmers in such configurations overall benefit in their propulsion (figures 3(b) and 4(b)), suggesting that a swarm of linear swimmers should propagate faster if the swimmers arrange behind each other or all side-by-side [43]. However, this result is sensitive to the relative swimmer positioning and orientation. In particular, if the swimmers are arranged on a line or side-by-side, but with opposing swimming directions, the swimmers will predominantly hinder each other's propulsion due to passive effects.
Comparing the active interactions using the map from force parameters to stroke parameters, we find that the results in both models generally differ, but agree for arbitrary swimmer positioning in the limit Γ = 6πηaω/k → 0, i.e. when the driving frequency in the FB model ω becomes small compared to the inverse viscous time t −1 V . In figure 5, we plot the total quadrupolar interaction along the swimmer axis (Δv par ) to second order in A and ξ, respectively, for two side-by-side swimmers in dependence of the phase shift γ, while varying the rescaled driving frequency Γ. It becomes apparent that for Γ approaching zero, the respective FB curve approaches smoothly the SB curve until they coincide in the limit Γ → 0. In contrast, for large values of Γ, the interaction effect in the FB model is out of phase by approximately π compared to the SB active interaction, which means that while active effects in one model boost the swimmer, they slow it down in the other model. This shows that for frequencies with Γ ≈ 1 or larger, the SB results are insufficient to describe the interaction of elastic swimmers and a theory accounting for the altered swimming stroke is inevitable. While figure 5 corresponds to the interaction of two pullers, qualitatively very similar results are found for the different other combinations of pushers and pullers. Figure 5 is hence representative for interacting linear swimmers, independently of the dipolar swimmer types.
To compare both models in figure 5, a certain stroke amplitude and phase between the two swimmer arms had to be assumed for each swimmer in the SB model. The driving forces for the FB swimmers were then adjusted for each choice of frequency such that at the single swimmer level the resulting stroke amplitude and phase would coincide with the before fixed stroke. We highlight that when comparing the interactions within two-swimmer systems in both approaches, the driving forces on each swimmer have to be adjusted such that they produce the SB model stroke when the swimmer is alone. If we, hypothetically, would adjust the driving forces such that the resulting swimmer strokes would coincide with the SB strokes while the swimmers interact, we would simply recover the SB system since fixed strokes together with the force-free and torque-free conditions already determine the whole system dynamics.
Analyzing the terms responsible for the active interactions allows to get an intuitive understanding for why they are different in both models. Active quadrupolar translation of swimmer I, to second order in A and ξ, results from the time-dependent oscillating dipolar flow field, scaling as A 1 or ξ 1 respectively, produced by swimmer II, interacting with the stroke of swimmer I. This flow field can be decomposed in a Taylor series around the position of swimmer I. The zeroth order, spatially constant term is associated with the same hydrodynamic force on all three beads, an effect averaging out over one stroke cycle due to the purely oscillatory nature of the flow field. At the next order in the expansion, the gradient of this instantaneous dipolar flow field, which swimmer I senses along its axis, is associated with hydrodynamic forces which expand and compress the swimmer arms and is responsible for the active interaction. In the FB model, the interplay of those forces with the springs and the viscous friction alters the swimming stroke by an additional term ∼ A 1 /r 3 . Since the swimmer's velocity is effectively quadratic in its stroke, cross terms of this additional term and the swimmer's original stroke ∼ A 1 r 0 yield the leading order active interactions with the correct scalings. In the SB model, the arm lengths are prescribed a priori such that the gradient of the time-dependent dipolar flow field evokes additional external forces necessary to maintain the prescribed stroke rather than an altered stroke. The active interactions then arise as a consequence of the additional forces [37].
From this, we can understand why both the SB and FB models in general differ in the active interactions, but become equivalent when in the FB model Γ = 6πηaω/k → 0. This limit can also be interpreted as the limit of k → ∞ when fixing ω and η, i.e. as the limit of very stiff springs. Consequently, to maintain a constant swimming stroke amplitude of a single swimmer, it is then necessary to scale up the driving force amplitude accordingly. A such FB swimmer with stiff arms will barely alter its stroke in response to compressing or expanding hydrodynamic forces produced by another nearby swimmer. Instead, the stiff springs will exert forces along the arms counterbalancing the hydrodynamic forces-and the mechanism of the active interactions becomes identical to the one for the SB swimmers.
Comparing to the literature, we find that the results (18) and (19) have not been reported yet. Other authors have obtained results for the quadrupolar active interaction scaling to leading order linear in L [37,44]. Again, we have reproduced such linear scaling in a calculation employing the altered perturbative scheme, which assumes the distance between the middle beads of both swimmers to be fixed to r throughout the swimming stroke. Conversely, we find in our calculation that the active quadrupolar interaction scales to leading order constant in L, which is in agreement with our numerical calculations (figures 2-4) and also backed up by the explanation of the active interaction in terms of time-dependent dipolar flow field, which yields a similar scaling. We are able to confirm that the second order quadrupolar translational interaction acts along the swimmer axis only, a finding which has been reported previously [37].
In contrast to the quadrupolar translation, the octopolar (∼ 1/r 4 ) translational interaction has in general both a component along the swimmer axis, which typically becomes important only at close swimmer separations since it decays quicker than quadrupolar effects, as well as a component acting orthogonal to the swimmer axis. As the quadrupolar translational interaction has no orthogonal component, a transition from the passive dipolar far field interaction to active octopolar interaction at intermediate separations can be observed for Δv ort in the side-by-side swimmer configuration (figures 4(d) and (e)) at a few tenths of swimmer lengths. For two equal swimmers and γ = 0, the active orthogonal octopolar component vanishes, therefore the transition is not observed in figures 4(c) and (f). In the SB calculation we obtain for the ξ 2 octopolar interaction of two side-by-side swimmers orthogonal to the swimmer axis Δv ort,SB Since all summands in (20) are dependent on δ, this term is purely active, consistent with the absence of even order ξ 2 flow fields and passive interaction for the linear swimmer geometry. This result is, to leading order in 1/L, in agreement with the corresponding formula from [39].
Comparing to the corresponding result for the FB interaction, Δv ort,FB we find that both results are indeed equivalent with respect to the mapping from the FB to the SB model. This means, this interaction effect is independent of the spring constant k in the FB model and equal to the interaction effect for two SB swimmers, assuming that the driving forces one each swimmer correspond to the same single swimmer stroke. This result holds for arbitrary swimmer configurations. In contrast, we find that the component of the octopolar translation along the swimmer axis is not equivalent in both models.

Active rotation
Active rotation is firstly observed at octopolar order (∼ 1/r 4 ) and order A 2 and ξ 2 . In the SB model, we find for two side-by-side swimmers where both the active and passive contributions are included. In the FB model, we again restrict the analytical result to two pullers in the same configuration, finding Employing the mapping from FB to SB parameters reveals that also this octopolar rotation is equivalent in both models. The active octopolar rotation scales to leading order as L 1 for L → ∞, consistent with the results reported in the literature [34,37,39]. In particular, we find that our result agrees algebraically in the leading order in 1/L with the result reported in [39]. Unlike the active contribution, the passive contribution to octopolar rotation scales to leading order as L 0 , as can be seen from (22) and (23). Therefore, the active component of the rotation is typically dominating over the passive one, and clockwise (CW) as well as CCW rotation can both be observed when varying γ or δ ( figure 4(b)) due to the sinusoidal dependence of active interaction on the phase shift.
Investigating the origin of the terms responsible for the active octopolar rotation at this order in the perturbative scheme, we find that three different terms contribute to it. First, swimmer I rotates with magnitude ∼ A 1 /r 4 in oscillatory fashion due to the curl associated to the time-dependent quadrupolar flow field produced by swimmer II. From the interplay of this time-dependent rotation with the swimmer's own self-propulsion, swimmer I experiences an active overall rotation ∼ A 2 /r 4 . Second, swimmer I moves within the instantaneous flow field produced by the other swimmer due to its A 1 r 0 original swimming stroke, experiencing the curl of this flow field at different positions. The inference of both the oscillating motion and the time-dependent curl of the flow field contributes to the active interactions as well. Third, swimmer I experiences hydrodynamic forces which would bend the swimmer, resulting from the instantaneous dipolar flow field produced by swimmer II. The counteracting forces exerted by the rigid joint, although themselves being torque-free, induce fluid flows, which, in interplay with the A 1 r 0 swimming stroke of swimmer I, give also rise to active rotation. All three of these mechanisms are independent of whether the swimmer arms are stiff or elastic, explaining why the active octopolar rotation is equivalent in both models.

Swimmer behavior for arbitrary positioning and long-term behavior
We now extend our considerations to two FB swimmers positioned arbitrarily in the two-dimensional space and the temporal evolution of their relative positions. For the latter, different types of long-term trajectories for two swimmers have already been classified in the SB model in dependence on the initial relative positioning [34,37]. Both authors investigating the long-term behavior have obtained broadly similar results for the long-term trajectories, although different driving protocols were employed. We extend these studies to FB swimmers with sinusoidal driving and show that the long term behavior is qualitatively the same as for SB swimmers [34,37] and is also robust with respect to varying the swimmer characteristics between puller-and pusher-type. This can be understood from the fact that the long-term behavior is predominantly a result of the self-propulsion and rotational interaction, which are to leading order equivalent in both approaches, and is only weakly influenced by translational interaction [34]. In contrast to previous works, we infer the long-term behavior from the instantaneous interactions evaluated in the parameter space of all possible relative swimmer configurations, allowing for a simple graphical representation of the system's dynamics.
The relative positioning of both swimmers is parametrized by their distance r, the orientation of the connection between the middle beads of both swimmers relative to the axis of swimmer I, ϕ, and the difference between the two swimmer orientations, θ II − θ I (figure 1). The equations governing the momentary evolution of a system of two swimmers, obtained by elementary geometry, are given by with ϕ = π − ϕ + θ II − θ I . Δv  denotes the corresponding time-averaged single swimmer speed. The first four terms on the right-hand side of (24) correspond to how ϕ changes due to the swimmer translation, while the last term accounts for the effect due to the rotation of swimmer I. To illustrate the behavior of both swimmers, we display the difference in the angular velocities of both swimmers, ΔΩ II − ΔΩ I , time-averaged over a complete stroke cycle in dependence of ϕ and θ II − θ I for a fixed distance r by the color in figures 6(a) and (b). The momentary time evolution of the system is illustrated by arrows corresponding to the vector field An arrow pointing to the left or right indicates that ϕ increases or decreases, and an arrow pointing upwards or downwards corresponds to θ II − θ I increasing or decreasing. When both ϕ and θ II − θ I change with time, the arrow is tilted such that the horizontal and the vertical component represent the corresponding rate of change. For better visibility the lengths of all vectors have been normalized. We obtain qualitatively similar plots for two interacting pushers or for combinations of a pusher and a puller. Figure 6 can therefore be considered as representative for the interaction of linear swimmers at intermediate distances up hundreds of swimmer lengths, independently of the dipolar characteristics. This effectively results from the fact that the linear swimmer with equal bead radii produces only a weak average dipolar flow field, as discussed before, and interaction at the distance of r = 30a considered here is dominated by quadrupolar effects. It is in contrast to swimmer models with dominant dipolar average flow fields, which often exhibit strong differences in the collective behavior of pullers and pushers [40,50,51].
The time evolution for two in-phase swimmers (γ = 0, figure 6(a)) and two swimmers out-of-phase by γ = π ( figure 6(b)) is dominated for a major part of the parameter space spanned by ϕ and (θ II − θ I ) by the self-propulsion terms in (24), which are linear in v 0,s . In these parts of the parameter space, we have |d/dt ϕ| |d/dt(θ II − θ I )|, associated with horizontal arrows. Exceptions from this are first the regions around the green dashed lines, corresponding to θ II − θ I = 2ϕ ± π, and second the region around the ϕ-axis, i.e. θ II − θ I = 0, corresponding to swimmers with parallel swimming direction. It is straight forward to verify that the first case corresponds to both swimmers being symmetric with respect to some axis in the x-y-plane, as sketched in the inset in the first panel of figure 6(c). In both of the two cases, the self-propulsion terms proportional to v 0,s in (24) cancel. Therefore, the swimmer behavior in these regions is dominated by rotational interactions. In the symmetric configuration (green dashed lines), in-phase pullers rotate towards each other (which can be seen from the vectors pointing down at ϕ = π/2, θ II − θ I = 0 in figure 6(a)), but rotate away from each other when out-of-phase, consistently with figure 4(b).
During swimmer interaction, not only the values of ϕ, θ I and θ II change, but also the distance r between both swimmers does. Varying the value of r, the angular velocities ΔΩ I and ΔΩ II scale as r −4 as long as r does not leave the regime dominated by the quadrupolar flow fields. Hence, the color structure in figure 6 stays constant while only the value corresponding to each color is rescaled. Conversely, the arrows align with the horizontal direction when r is increased, since the horizontal component (24) scales to leading order as r −1 due to v 0,s being independent of r. In contrast, the vertical component (25) scales as r −4 .
Although figures 6(a) and (b) correspond to constant r and cannot account for the swimmers coming closer or separating, we are able to reproduce the different types of long-term trajectories (figure 6(c)), which have been similarly reported for two interacting linear SB swimmers with initially parallel swimming direction (figures 4 in [34] and figure 3 in [37]). In figure 6(c), the middle bead of each swimmer is taken as reference for the swimmer position. Since the translational interaction is small compared to self-propulsion, both swimmers are in general oriented tangential to their trajectories, swimming from left to right.
We first consider two swimmers with zero phase shift γ ( figure 6(a)) and classify the different types of long-term trajectories observed when the swimmers have initially parallel swimming direction (θ II − θ I = 0) but are positioned arbitrarily otherwise, corresponding to varying the initial value of ϕ. Two approximately collinear swimmers, i.e. ϕ ≈ 0, swim on oscillatory (label O in figure 6) trajectories. Mapping this behavior to the phase space of ϕ and θ II − θ I , it corresponds to the curl around ϕ = 0, θ II − θ I = 0 found in figure 6(a). With increasing ϕ, we next enter the repulsive (label R) regime in which the swimmers rotate away from each other and thus separate such that rotational interaction dies out and the swimmers self-propel on straight trajectories away from each other. Increasing ϕ further, still keeping the swimmers initially parallel, the swimmers typically move on parallel trajectories (label P), associated to the curl close to ϕ ≈ π/3, θ II − θ I = 0. Since this curl is associated to positive values of ϕ only the swimmer trajectories do not cross, distinguishing it from the oscillatory case. For ϕ close to π/2, the swimmers rotate towards each other, effectively attracting each other (label A), until the beads come so close that the Oseen approximation is not valid anymore. This behavior corresponds to the area around the green dashed lines in figure 6(a).
Similarly, also the types of trajectories of two out-of-phase swimmers are accounted for in the respective plot of the ϕ-(θ II − θ I ) space ( figure 6(b)). Besides the types of trajectories already mentioned, we also observe trajectories where the swimmers first rotate towards each other, their trajectories cross and the swimmers subsequently repel (label O R ). Also for this behavior, the qualitative trajectory in terms of ϕ and θ II − θ I can be depicted in figure 6(b).

Discussion and conclusions
We have calculated the scaling function, the strength and the direction of the hydrodynamic interactions between two linear three-bead swimmers. Both the FB and SB models were considered perturbatively, and the findings were very favourably compared to direct numerical results. While this framework is applicable to arbitrary positioning of the two interacting swimmers, the most simple cases of collinear and side-by-side swimmers have been used to illustrate the main results of the paper. Due to the simplicity of the swimmer design, these results may be relevant in a number of systems and find application well beyond the physical context presented herein.
At the leading quadratic order in the swimmer actuation, the time-averaged flow field produced by a linear three-sphere swimmer, defining passive hydrodynamic interactions, is proportional to its velocity (section 3.1 and appendix D). The proportionality coefficient depends solely on the swimmer geometry. This novel result, valid for FB as well as SB swimmers, allows us to infer the passive interactions to second order in the swimmer actuation. The latter are given by the time-averaged flow field produced by all nearby swimmers, and can be directly inferred from the geometry of the system. The first consequence of this finding is reflected in the identical passive interactions of FB and SB swimmers under the constraint that the individual devices display identical actuation. Another consequence is the enhancement of the swimming velocity of both swimmers placed one behind another (figure 3), or side-by-side (figure 4), compared to the velocity of a single swimmer. This cooperative effect for pairs of swimmers is independent of the details of the driving and is found similarly for pullers and pushers. It relies on the passive quadrupolar interactions, which dominate translations in a major part of the parameter space.
Active hydrodynamic interactions, which are the result of the interference of the time-dependent swimming strokes of both swimmers, are particularly important for rotations. Nonetheless, active torques ( (22) and (23)) as well as translational interactions orthogonal to the swimmer axis ( (20) and (21)) are equivalent to the leading octopolar order in the SB and the FB model. Differences are, however, expected in the higher orders. Notably, the active translations along the swimmer axis differ already in the leading quadropolar order ( figure 5). This shows that the interaction is in general not conserved with respect to the mapping between FB and SB swimmers and that the deformation of one swimmer due to the presence of others significantly alters the change in propulsion speed that a swimmer experiences, depending on the design. Nonetheless, the FB and SB models become equivalent when the driving frequency is small as compared to the inverse viscous time, which is meaningful only in the FB model. Consequently, the SB interaction should be regarded as the special case of the FB interaction in the limit of small driving frequencies or high swimmer stiffness.
We, furthermore, find that the long-term behaviors for FB and SB swimmers agree qualitatively, which is a result of the equivalence of the octopolar rotational interactions in both models. In the long-time limit, the planar two-swimmer systems are, furthermore, independent of the pusher/puller characteristics of either swimmers up to separations of hundreds of swimmer lengths. Actually, the long-term trajectories can be qualitatively inferred directly from the instantaneous behavior of both swimmers, as evaluated numerically in the space of all relative swimmer configurations ( figure 6). While these considerations here have been restricted to planar swimmer configurations, the framework permits arbitrary orientations of the swimmers. In this case, we expect no qualitatively new behavior in the instantaneous swimmer interactions, however, three-dimensional initial configurations might indeed give rise to interesting new phenomena for the long-term behavior. Preliminary numerical calculations suggest for instance that the oscillatory regime shown in figure 6 will transition into a helical one, when one of the swimmers is tilted out of the common plane of both swimmers. A comprehensive analysis of the three-dimensional trajectories is however beyond the scope of this paper and will be addressed in a future work.
We have demonstrated that active interactions are of great relevance for the three-bead swimmer, which can be seen as a toy model for swimmers which self-propel by periodic shape changes. Therefore, active interactions might also be of importance for biological microswimmers which exploit shape changes for their locomotion. One of the most prominent examples for such a species is the flagellated algae Chlamydomonas reinhardtii, which moves using a breast stroke-like beating pattern [52]. By adapting the swimmer geometry, the framework presented here can be applied with ease to a triangular three-bead swimmer, a geometry that has already been successfully used to model Chlamydomonas [53][54][55]. The role of the active interactions for the collective behavior of many such algae cells is therefore an interesting question which now can be addressed explicitly.
Another interesting question is the role of the elasticity in biological microswimmers. Experimental studies have shown that algae cells such as Chlamydomonas [56], but also ciliated microorganisms [57] do indeed change their stroke in dependence of the environment, suggesting that a purely SB approach might fail to fully capture features as the swimmer behavior in swarms or in spatially varying fluid properties such as viscosity. We showed that FB swimmers can adapt their stroke to neighboring swimmers, and that the strength of this adaptation can be manipulated by the swimmer elasticity. As shown in this paper, the degree of adaption is an important parameter controlling the active interactions. However, it is yet unclear if such a simple mechanistic approach can capture all relevant facets of the swimmer dynamics observed experimentally. An appropriate description might comprise both FB and SB elements, or even more complex feedback mechanisms. However, the fact that swimmer interactions do depend on the swimmer elasticity points towards the possibility of employing bead-spring swimmers as a simplistic toy model in order to gain new insights in the mechanisms underlying self-locomotion of groups of biological organisms.
Indeed, this work establishes a toolbox for a better understanding of the hydrodynamic interactions in swarms of microswimmers, which can be studied both numerically and analytically. Naturally, increasing the number of swimmers as well as employing a more precise hydrodynamics will be accompanied by growing computational costs. However, while it still accounts for the full hydrodynamics at a tunable level of precision, our framework should be still computationally less expensive than a simulation method which resolves the hydrodynamics at a microscopic level, such as the lattice Boltzmann method. This should be particularly notable in dilute swarms where the restriction to pairwise microswimmer interactions in determining the swarm behavior seems natural. In this case, the analytical results for two swimmers presented in this work can be used directly, appreciably speeding up the calculations, while maintaining accuracy.
The theory developed herein should be also applicable in dense swarms, where hydrodynamic interactions were also shown to be important [58]. Even in these conditions, it may be sufficient to restrict the analysis to two-swimmer interactions. Namely, three-body effects decay at least with the fifth inverse power of the typical swimmer separation and should thus become important only when the swimmers approach each other closer than approximately three to four swimmer lengths. At these densities, the numerical approach presented here can be employed. This will also require a higher order hydrodynamic theory, as the application of the Rotne-Prager tensor. This could unfortunately increase the computational costs of the presented framework, making the calculations similarly demanding as in mesoscopic methods. In this dense regime, however, based on the observation that swimmers in pairs typically speed up when they are arranged collinearly or side-by-side, we predict that small swarms will be able to self-propel faster than a single individual can. Larger swarms with in-phase swimmers arranged on a rectangular lattice have shown that such configurations are prone to instability and collisions between the swimmers [28]. This result is rooted in the tendency of such swimmers to rotate towards each other and accelerate. Yet, it is an interesting question whether stable or at least metastable swarming states exist or if a more advanced feedback mechanism is required, a problem which will be addressed in a future work.
All dashed variables denote rescaled quantities with respect to the bead radius a and the viscous time t V [9], τ := t/t V denotes the rescaled time t and := A/(ka). In particular, we have introduced (n · d)-dimensional bold vectors for the rescaled position and the external driving forces 3) The driving forces as defined in (4) have to be treated in the perturbative formalism as dependent on the bead positions because they act along the swimmer orientation vectors n s . Denoting the total spring force on each bead sj by as the vector of all rescaled spring forces on all the beads. The mobility matrix μ (R ) with (n · d) × (n · d) components is defined by n × n blocks with in the diagonal blocks 1/(6πηa) the self-mobilities and in all other blocks the corresponding Oseen tensor (see (8) in [9]). In addition to the driving forces, which act along the respective swimmer axis in the form specified in (4), we include anti-bending forces for each swimmer and at each order in . They are denoted, in rescaled form, by E (i) A (R , τ ) for forces at order i . The form of the anti-bending forces on each bead of a swimmer is derived from a potential defined to be the angle enclosed between both swimmer arms, In contrast to the potential φ s defined in appendix A, which attains a minimum at κ s = 0 corresponding to zero anti-bending forces, the potential Φ s is associated with non-zero forces that would bend the swimmer s even if it has linear shape. We then define the anti-bending forces as with A bend,(i)

II
(τ ) a priori undetermined time-dependent prefactors. This definition ensures that the overall force and torque exerted by the anti-bending forces on each swimmer vanish.
At order 1 we typically observe bead oscillations with the base frequency [9], hence it suffices to assume for the explicit time dependence The prefactors will be carried along in the calculation of the swimmer behavior at each order i and then determined a posteriori such that they effectively cancel out any bending of the swimmers. To solve (B.1) by perturbation theory around = 0 and q = 0, corresponding to zero driving forces and infinitely separated swimmers, we introduce the displacement of all beads with respect to their initial position, ξ (t) := R (t) − R eq , following the notation in [9], with R eq the vector of all rescaled initial bead positions. We assume that in this configurations both swimmers are in mechanical equilibrium, i.e. the distance between neighboring beads is L. We then expand the displacement as a power series in both and q, The first index of ξ '(i,j) is associated to powers of and the second index is associated to powers of q. For simplicity, we subsequently omit the time-dependence of ξ and R . We proceed by expanding μ (R ), G (R ) as well as the driving and anti-bending forces in a Taylor series around the equilibrium configuration of both swimmers, R eq , and insert the expansion (B.10) for the bead displacement (see [9] for details of the expansion). At each order of i q j , the equation of motion can then be cast into the form d dτ Here, K (0,0) denotes the q 0 component of the -independent K defined by K := μ (R eq ) · ∇ R G (R eq ), (B.12) with ∇ R the gradient with respect to all (n · d) bead positions, and S (i,j) (τ ) a source term which pools all terms independent of ξ (i,j) . Since the expansion of the bead displacements (B.10) contains only terms with positive exponents in and q, the source term S (i,j) (τ ) can only contain components of the displacement ξ (i ,j ) with i i, j j, but not i = i and j = j. The source term S (i,j) (τ ) composes at each order i q j of a contribution resulting from the expansion in and a contribution resulting from the expansion in q, q (τ ). The first contribution is given similarly as in the case of a single swimmer (equations (21) and (22) in [9]) for the first two orders in as, Here we use the notation S (i,all) = j q j S (i,j) . The expressions for S (i,j) (τ ) are then obtained by expanding the above expressions by powers of q. The second contribution S (i,j) q , due to the application of perturbation theory in q, is given at each order i q j by Since S (i,j) (τ ) depends only on the components of the displacement ξ (i ,j ) with i i, j j but not i = i and j = j, an ascending iteration scheme through the orders in and q allows to directly calculate at each order first the source term and second from this the bead displacement by solving (B.11). Being interested in the swimmer behavior up to e.g. order 2 q 4 , a suitable iteration scheme is given by To solve (B.11) at each order, we decompose the respective source term into an oscillating component and a constant component, with f i an integer and S (i,j),sin f , S (i,j),cos f , S (i,j),const prefactors to be determined by decomposition. Due to the linearity of (B.11), one can calculate the solution for each component separately and superimpose those solutions to obtain the final result for the bead displacement. For the oscillating source terms, the corresponding solution is given by [3,9] ξ (i,j) = Here, 1 denotes the unit matrix on the (n · d)-dimensional space of all bead positions. The constant contribution to the source term gives rise to both the translational and rotational velocities of both swimmers as well as to deformations. Solving (B.11) becomes simple by decomposing the source term in the eigenbasis of K (0,0) . By its definition, K (0,0) is independent of q and thus does not couple between both swimmers. As an example, for θ I = θ II = 0 the matrix K (0,0) is given by Similarly to the case of a single swimmer [9], the components of the constant source term associated with eigenvalue zero in the eigensystem of K (0,0) correspond to translation or rotation of one of the swimmers. In contrast, the components associated with negative eigenvalues correspond to swimmer deformation along its axis, i.e. the average arm lengths of the swimmer differ from the spring length in the mechanical equilibrium. Positive eigenvalues are excluded due to the stability of the mechanical equilibrium of the swimmer [9].
In the eigenbasis of K (0,0) , (B.11) assumes for the constant source term the form with λ (δ) , δ = 1, . . . , n · d the eigenvalues of K (0,0) . The overline denotes expression of a vector with respect to the eigenbasis of K (0,0) and the additional upper index δ the respective component of the vector. The translational and rotational velocities of each swimmer are thus given directly by the components of the constant source term corresponding to zero eigenvalues, since in this case λ (δ) = 0. In order to determine these components, we use a change of basis from the eigensystem of K (0,0) to standard coordinates, M , given by With Q = diag(1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0) the diagonal matrix with unity in the first eight entries and zeros everywhere else, translations and rotations of both swimmers are obtained as A suitable decomposition of these bead velocities into translation along the swimmer axis, orthogonal to it and rotation of each swimmer yields the respective velocities. Taylor expanding these velocities by powers of q then allows to identify quadrupolar, octopolar etc interaction effects.
A posteriori to the calculation of the swimmer behavior, we determine the amplitudes of the anti-bending forces at each order i by imposing with ξ (i,all) = j q j ξ (i,j) . Performing the steps described in this appendix involves tedious algebraic operations, which is why we use the software Mathematica [45] to perform the actual calculations.
Appendix C. Details on the stroke-based perturbative calculation C.1. Flow field of a single stroke-based linear three-sphere swimmer To calculate the flow field produced by a single linear three-sphere swimmer, we firstly calculate the trajectory of each bead and from this the velocity field in the fluid. We hereby follow the calculations performed in [39,48].
Similarly to the case of two interacting swimmers, the relation of the bead velocities and the forces necessary to enforce the prescribed swimming stroke is given by (12). In the case of one swimmer aligned with the y-axis, a one-dimensional framework with R(t) = (y 1 (t), y 2 (t), y 3 (t)) and is sufficient to calculate the bead trajectories. In contrast to the case of two interacting swimmers, here the mobility matrix μ(R) is determined already by the prescribed stroke, simplifying the calculation. With the two constraints resulting from the prescribed arm lengths, as well as the force-free condition, F 1 (t) + F 2 (t) + F 3 (t) = 0, the system of equations (12) is closed and can be solved using Mathematica [45]. Here, the dot denotes the derivative with respect to time. With the prescribed stroke inserted, it is straight-forward to obtain the bead positions from the respective velocities by integrating with respect to time. The flow field produced by the swimmer, u( r, t), is then given by [39] u( r, t) with x i (t) = (0, y i (t)) and F i (t) = (0, F i (t)) the extensions of scalar quantities to two-dimensional vectorial quantities. Time-averaging and expanding the flow field with respect to r allows to decompose the flow field obtained into dipolar (∼ 1/| r| 2 ), quadrupolar (∼ 1/| r| 3 ), . . . components.

C.2. Interaction of two stroke-based linear three-sphere swimmers
In this section, we elaborate on the details of the perturbative calculation for two interacting swimmer in the SB model, which is outlined in section 2.2. A major complication in comparison to the calculation of the behavior of a single swimmer is that here the mobility matrix μ(R) does not only depend on the arm lengths within each swimmer, which are prescribed, but also on the orientation of each swimmer and the relative positioning of both swimmers. Since the latter degrees of freedom are a priori unknown, we must use a perturbative approach in the inverse swimmer separation q to calculate their behavior. To speed up the calculation, we additionally perform an expansion also in the actuation amplitude ξ.
With the constraints on the swimmers' linear shape as well as with the prescribed arm lengths, the positions of all beads are parametrized by the six remaining degrees of freedom by R I2 , R II2 , θ I and θ II . We expand each of those quantities in a power series with respect to ξ and q: where the first upper index in brackets corresponds to powers of ξ and the second index to powers of q. Note that we omit the time-dependence of the bead positions, swimmer orientations and swimmer arm lengths for the sake of brevity. R init I2 , R init II2 , θ init I and θ init II denote the initial conditions for the two swimmers at t = 0. The positions of the remaining beads are then given for each swimmer s by With this expressions at hand, we are able to exploit (12) in order to solve for the forces F. In order to keep the calculation simple, we expand (12) with respect to q and find F (all,0) = μ (all,0) d dt R (all,0) , F (all,1) = μ (all,0) d dt R (all,1) − μ (all,1) F (all,0) , ( C . 5 ) F (all,2) = μ (all,0) d dt R (all,2) − μ (all,1) F (all,1) − μ (all,2) F (all,0) , where the first index 'all' indicates that the corresponding quantity is not expanded with respect to ξ. A subsequent expansion of the right-hand sides of (C.5) allows to find expressions for the forces expanded with respect to both q and ξ. The six remaining degrees of freedom, R init I2 , R init II2 , θ init I and θ init II , are then determined by the force-free and torque-free conditions (10) and (11). Expanding these equations by powers of ξ and q, we notice that the force-free condition at each order ξ i q j involves only F (i,j) whereas the torque-free condition involves also components F (i ,j ) with i i and j j. Therefore, we employ again an ascending scheme through the orders in ξ and q given in this case by q 0 ξ 1 → q 0 ξ 2 → q 1 ξ 1 → q 1 ξ 2 → q 2 ξ 1 → q 2 ξ 2 → q 3 ξ 1 → ... Using this ascending scheme, all lower order components of the forces or R I2 , R II2 , θ I , θ II required at order ξ i q j have already been computed. Consequently, the force-free and torque-free condition at each order ξ i q j , with all lower order components of F and R I2 , R II2 , θ I , θ II substituted by their explicit result, become algebraic equations for˙ R II . By solving these equations and time-integrating the results found we obtain the bead displacement at each order. A suitable decomposition into translational velocity along the swimmer axis, orthogonal to it and into rotational velocity yields the results for the different components of the interactions. Due to the increasing number of terms, we use Mathematica [45] to perform the calculation.

Appendix D. Proportionality of ξ 2 average flow field and swimming velocity of a single swimmer
We show that the amplitude of each ξ 2 component (dipolar, quadrupolar, . . . ) to the time-averaged flow field produced by the linear swimmer can be expressed as a purely geometric prefactor times dt[ξ 1 (t)ξ 2 (t) − ξ 1 (t)ξ 2 (t)], which represents the area in the swimmer's configuration space enclosed by its trajectory.
As noted by Golestanian and Ajdari [48], the instantaneous velocity of each bead of a single swimmer along its axis, v i , is linear in the velocities of the arm extensions and thus can be expressed as with b, c indices iterating over the swimmer arms and V, W general geometrical prefactors. Repeated indices are summed over and the dot denotes derivation with respect to time. We assume the swimmer is oriented along the y-axis. Due to the linearity resulting from the Stokes equations, the fluid velocity u( r, t) at some arbitrary position r in the absolute coordinate system and time t can then be written as a sum over all beads, u( r, t) = iÛ ( r − R i (t)) · (0, v i (t)), (D.2) with (0, v i (t)) the 2D velocity vector of bead i andÛ some tensor which depends on the vector connecting the position considered to the position of bead i. From integration of (D.1) we know that up to order ξ 1 , i.e. linear in either ξ b (t) orξ b (t), the position of each bead i is a constant plus a term linear in ξ b (t) with b = 1, 2. Using a Taylor expansion around r, thus alsoÛ( r − R i ) can be written, up to order ξ 1 , as a constant plus a term linear in ξ b . Hence, combining (D.1) and (D.2), we find that all ξ 2 contributions to u( r) can then be written as u( r, t) = C bc ( r) ·ξ b ξ c , ( D . 3 ) with C bc ( r) a prefactor which, besides the dependence on the position r in the absolute coordinate system, depends only on the geometry of the system. The only combination ofξ and ξ which does not vanish when averaged over one stroke cycle isξ 1 ξ 2 − ξ 1ξ2 . This proofs our claim. We point out that in contrast to the velocity of each bead, the flow field at a fixed position in the absolute coordinate system does not depend linearly onξ b (t) when we consider terms of order higher than 2 in ξ. Since the swimmer self-propagates in the absolute coordinates and therefore relative to r, the vector connecting r to R i will at order ξ 2 also contain terms of the form t 0 [ξ 1 (t )ξ 2 (t ) − ξ 1 (t )ξ 2 (t )]dt . Therefore, the instantaneous flow field at order ξ 3 will also contain terms of the form t 0 [ξ 1 (t )ξ 2 (t ) − ξ 1 (t )ξ 2 (t )]dt ξ b (t) and the expansion (D.1) is not valid for components of the instantaneous flow field.