Understanding the Onset of Oscillatory Swimming in Microchannels

Self-propelled colloids (swimmers) in confining geometries follow trajectories determined by hydrodynamic interactions with the bounding surfaces. However, typically these interactions are ignored or truncated to lowest order. We demonstrate that higher-order hydrodynamic moments cause rod-like swimmers to follow oscillatory trajectories in quiescent fluid between two parallel plates, using a combination of lattice-Boltzmann simulations and far-field calculations. This behavior occurs even far from the confining walls and does not require lubrication results. We show that a swimmer's hydrodynamic quadrupole moment is crucial to the onset of the oscillatory trajectories. This insight allows us to develop a simple model for the dynamics near the channel center based on these higher hydrodynamic moments, and suggests opportunities for trajectory-based experimental characterization of swimmers' hydrodynamic properties.

The locomotion of self-propelled particles (swimmers) typically occurs at boundaries or under confinement. Accurately describing the effect of confinement on swimmers is therefore of significant interest to understanding the behavior of microorganisms and artificial swimmers. In modelling these systems, hydrodynamic interactions (HIs) are often ignored, which is a valid approximation in some cases, such as when microbial swimmers' run-andtumble dynamics dominate [1]. However, HIs can play an important role, e.g., see Refs. [2][3][4][5][6], and therefore cannot be a priori ignored in modelling. Recent experiments on self-phoretic colloidal swimmers have shown that their orientation is strongly influenced by HIs due to the presence of a wall [7]. However, there is ongoing debate on the importance of near-wall effects and the level at which to truncate the hydrodynamic moment expansion [8].
A specific example of this are the helical and oscillatory trajectories of single swimmers in confining geometries as observed experimentally by Jana et al. [5] and in simulations [9][10][11]. Such oscillatory trajectories appear to be common place, having been reproduced by many models, and independent of specific swimmer type. However, a physical understanding of these oscillations remains wanting. It is indisputable that the oscillations do not arise simply from the lowest order hydrodynamics, which result in direct attraction to surfaces [12], while the inclusion of higher-order modes can lead to more complex behavior [8,13]. Although observed in confined quiescent fluids, these oscillatory trajectories are reminiscent of those observed in the rheotaxis of swimmers subjected to external flows [14], which result primarily from the interplay between the flow and persistent particle motion due to self-propulsion. Zöttl and Stark indicate that near-field lubrication theory can be used when there is no externally applied flow to describe such trajec- * jgraaf@icp.uni-stuttgart.de tories [15,16]. Yet, the observations of Zhu et al. demonstrate that oscillatory trajectories arise in a channel that is three times as wide as the self-propelled particle. Additionally, the trajectories of squirmers close to single walls in quiescent fluid show oscillations [17], which have been explained by the competition between far-field HIs and short-ranged wall-swimmer potentials [18]. Thus, there is a clear need to establish to what extent the observation of oscillatory trajectories in systems with confinement originate from a near-or far-field effect and, in conjunction, to assess the importance of higher-order hydrodynamic modes.
In this manuscript, we demonstrate that the onset of time-varying oscillatory trajectories in systems confined within a channel and without external flow can be wellunderstood using far-field theory [19]. We investigate the specific case of two parallel infinite plates that enclose the fluid and a single rod-shaped swimmer, using our lattice-Boltzmann (LB) 'raspberry' force/counter-force formalism [20] (Fig. 1; insets). We have previously shown that the rod-shaped LB swimmers have well-defined higherorder hydrodynamic moments [20]; see Table I for representations of the first five moments. These simulations conclusively show that a puller-type rod that starts far from the wall but off-center follows a damped oscillatory trajectory towards the middle of the channel, whereas a pusher-type rod moves between the walls along a sinusoidal path with increasing amplitude. Surprisingly, the oscillations are observable even for plate separations as great as ten times the length of the rod. We explain these observations within the framework of our far-field hydrodynamic theory: the dipole and octupole moments induce hydrodynamic forces towards the center (puller) or towards the walls (pusher), while the quadrupole moment causes pure oscillatory motion. The oscillatory trajectories within plate confinement thus provides an indirect means to characterize the hydrodynamic properties of swimmers, which would grant access to moments beyond those that can be obtained from lattice swimming [21] or tracer paths [20]. We consider two raspberry swimmers (rod and cylinder) in the main text to study the movement of shape-anisotropic swimmers under confinement using our ESPResSo LB implementation [20,22]. Their construction and characterization in terms of hydrodynamic moments, as well as the fluid and coupling parameters, are introduced in Refs. [20,23] and detailed in Appendix A 1. Our swimming model's essential aspect is that a force is applied to the body, consisting of many fluid-particle coupling points, and the system is made force free by applying an equal and opposite force to the fluid, see the insets in Fig. 1. This coupling gives rise to a series of hydrodynamic modes for anisotropic particles [20].
These raspberry particles are placed in an LB fluid between two parallel (no-slip) bounce-back plates, with normals in theẑ direction, separated by a distance H. The fluid domain is periodic in the other two (xy) directions. The vertical position of the swimmer's center of mass (CM) is indicated using z ∈ [−H/2, H/2], with z = 0 the middle of the channel. Lateral displacement is given by x and measured from the swimmer's initial position (x = 0) -our trajectories are straight in the xy-plane. Finally, the angle of the swimmer's directorp (which points along the main axis) with the plate normalẑ is given by φ ∈ [−π/2, π/2], with φ = 0 swimming parallel to the plates. To prevent the swimmers from penetrating the wall, we imposed a short-ranged (almost hard) Weeks-Chandler-Anderson (WCA) interaction between the raspberry swimmers and the walls (Appendix A 2). This wall-swimmer interaction is necessary as our LB algorithm does not explicitly account for near-wall lubrication corrections [23]. All of the results shown in the main text employ a WCA diameter d = σ, with σ the LB lattice size. We limit the swimming speed to ensure the low Reynolds number regime, Re < 0.01. Figure 1 and the supplemental movies (not included in arxiv version) demonstrate the onset of oscillatory trajectories. These are representative sample swimmer trajectories, where the swimmers start off-center and oriented parallel to the plates. Both the rod and cylinder models of pushers and pullers display time-varying oscillatory behavior. In the specific case of our rods, the wavelength of the oscillations is λ ≈ 4H. All pushers move towards the wall and the pullers move towards the center of the plates. After only a few periods, these pullers move along the centreline of the channel and these pushers have arrived in the near-wall region, where swimmer specific details and lubrication corrections would be required to accurately predict dynamics. Oscillations are observed for all cylinder and rod swimmers in plate separations that we could simulate (5σ ≤ H ≤ 50σ). The rod is ∼ 5σ in length, thus the oscillatory trajectories arise in systems with a channel height to particle size ratio up to at least 10.
To verify the generality of the initial oscillations, we considered several initial positions z and orientations φ for rod pusher and puller swimmers. We found that depending on the type of swimmer and its initial position/orientation, several oscillations in the physical regime can be observed, before near-wall effects cannot be ignored. We further showed that oscillations for rodlike swimmers appear for a large range in rod aspect ratios (Appendix A 6). At long times the LB pusher rods display a limit cycle, whereas the pusher cylinder does not. To what extent such a limit cycle ( Fig. 1a; solid red curve) or sliding dynamics ( Fig. 1a; dashed blue curve) might be physical is not considered here, as algorithmic artifacts close to the walls impact the near-wall dynamics. Appendix A 5 provides a discussion of these limitations and this work does not confirm their physical nature [11]. However, our results establish that the initial oscillations before the rod comes close to the wall (a proximity of ∼ 2σ) are physical. It is this onset of oscillatory trajectories that we concern ourselves with here and subject to theoretical analysis in the following.
We model the raspberry swimmers theoretically as el-lipsoids with aspect ratio γ, position r and orientation p = (cos φ, 0, sin φ). Due to its motion, the swimmer generates a flow field u, which we define in terms of a multipole expansion of the Stokeslet flow solution. Spagnolie et al. argue that far-field HIs give surprisingly accurate results, when compared to theory that includes a finite-size correction to more accurately account for nearfield effects, even for small swimmer-wall separations [8]. Hence, the flow at position x generated by the force-free and torque-free swimmer is Here, u D is the Stokes dipole that models the force balance between propulsion and drag, u Q is the quadrupole that represents the fore-aft asymmetry of the propulsion mechanism, u SD is the quadrupolar source doublet that is associated with the finite size of the swimmer, and u O1 and u O2 are the two octupolar terms (Appendix A 1). The shape of these moments in bulk is shown in Table I. Note that this is a point-based expansion, which should not be confused with the squirmer expansion for finitesized spheres; in the far-field these expansions can be mapped onto each other. The effect of the confining walls (two parallel noslip plates) is now accounted for by the method of images, where we truncate the approximation after four image systems on each side of the microchannel. Subsequently, the flow wall-induced flow advects and reorients the swimmer according to the Faxén relations, resulting in the translational and angular velocities v HI and Ω HI (Appendix B 1). In Ref. [19] details are given of the procedure by which to obtain these velocities in terms of the multipole coefficients. The swimmer's equations of motion are given bẏ where v s is the autonomous swimming speed, and the velocities v HI and Ω HI are functions of the multipole coefficients.
To predict the swimmer dynamics of Fig. 1 theoretically, we integrate the equations of motion (2). Here, we use the same swimming speed and initial conditions as in the LB simulations, but we allow the multipole coefficients to vary about their measured values. We can thus fit the multipole coefficients via a least-squares method. To obtain the best agreement with the measured trajectory, we used the four initial oscillations (Appendix C 1). The 3 rd and 5 th columns of Table I show the multipole coefficients found in this manner for swimmers of the rod and cylinder type, respectively. Using only a single oscillation leads to a change in the fitted values of ∼ 20%, showing our method to be robust and requiring only fragments of a trajectory to be effective. In addition, we verified that the result of the fitting was independent of the height H of the channel, eliminating the possibility of boundary artifacts. In our previous work [20], we  I. Multipole moments of the swimmer-generated flow field for the two swimmer types: the rod and the cylinder. The columns labelled 'LB' provide the values measured in our previous study by means of Legendre-Fourier decomposition in a close-to-bulk system with periodic boundary conditions [20]. The columns labelled 'theory' provide the moments fitted from the trajectory in our confining geometry by using the theoretical model (2). Values are given in LB simulation units, and the positive/negative signs correspond to pusher and puller swimmers, respectively. The bottom row shows representations of the flow field of the first five hydrodynamic moments in bulk: dipole κ, quadrupole ν, source dipole µ, source octupole o1, and octupole o2. The arrows are stream lines and the colors indicate flow away from (red) or towards (blue) the swimmer.
obtained the multipole coefficients directly from the flow field of the swimmers in our LB simulations by means of projection via a Legendre-Fourier decomposition. These values are listed in the 2 nd and 4 th columns of Table I, respectively. The projection was carried out in the absence of confinement, using a large simulation box with periodic boundary conditions, for which the finite-size effects differ strongly from those of the confining geometry. There is excellent correspondence between the two measurements of the hydrodynamic moments for both swimmer shapes. This demonstrates the applicability of far-field theory to describe the onset of the observed oscillatory trajectories. The far-field result is accurate until the swimmer-wall distance becomes too small. Let us now focus on the general features of the theoretical model and analyze the impact of the various hydrodynamic moments on the motion of the swimmer. Firstly, our calculations confirm that the pusher swimmer (κ > 0) undergoes oscillatory trajectories that move away from the center of the channel, and pullers (κ < 0) converge towards the centerline. However, oscillations about the center only occur if the quadrupolar terms are included, and the oscillation wavelength decreases with the associated quadrupolar coefficients ν and µ. A spherical swimmer with ν = µ = 0 [20] does not display such oscillations. The octupolar contributions further control the damping and growth of the trajectories, where the positive signs of o 1 and o 2 correspond to motion towards the boundaries. The aspect ratio γ leads only to a secondorder correction in the theory. That is, the hydrodynamic moments dominate the dynamics of the swimmer, therefore rods with different aspect ratios still show similar oscillations (Appendix A 6).
The dynamical system can be understood further by considering the motion of the swimmer in phase space. Due to the translational invariance in the x and y coordinates, the equations of motion can be reduced to two coupled first-order PDEs in (φ,z) space, next to the uncoupled equation for the x coordinate. Figure 2 shows the LB swimmer trajectories in phase space, superimposed with the theoretical model, where the fitted multipole moments in Table I have been used. The dipolar term leads to a star-type fixed point (curves radiating from a point) at the origin, that is stable for pullers and unstable for pushers. The oscillatory motion due to the quadrupolar contributions corresponds to a circle-type phase-space trajectory (closed loops around a point) centered on the origin. Together the dipole and quadrupole produce a spiral. For pushers, the trajectories spiral outwards (Fig. 2a), and inwards for pullers (Fig. 2b). The theoretical predictions do not show a limit cycle in Fig. 2. Both the far-field framework and the LB method are unable to adequately capture hydrodynamic interactions in the near-wall region and further study of this region, where both lubrication corrections [15] and short-ranged potentials [18] can play a role, is required.
Our result shows that movement of a swimmer under confinement can in principle be used to quantitatively determine the hydrodynamic moments, even up to the octupolar moment as shown here. Specifically, about one period is the minimum path length required to fit these modes to within 20%. This suggests that our method has applicability to experimental systems where thermal noise and tumbling can effect the trajectory. The presence of these sources of noise would require ensemble averaging trajectories in (φ, z) space, which can then be fitted using our procedure. Noise also implies that parts of the trajectory will occur many times during measurements, meaning the near-wall dynamics in which we observed a limit cycle, does not play an important role. One simply averages many different trajectories away from the wall to improve the fitting statistics.
A physical intuition for the onset of oscillatory swimming can be distilled from the LB simulations and farfield hydrodynamic description by considering the trajectory of a swimmer initially set at z 0 near the centerline and oriented parallel to the walls (Fig. 3). Since our raspberry swimmers have large quadrupolar moments (ν ∼ 10 −1 ), we first consider only the flow fields generated by the positive quadrupole. This flow serves to rotate the swimmer away from the nearest wall, as shown schematically in Fig. 3a. The continual rotation away from the nearest wall establishes the oscillations. By linearly expanding the equations of motion (2) about the centerline, the micro-swimmer dynamics can be captured by a linear system of coupled differential equations (Appendix C 2). Whenever there is only a quadrupolar flow field, the trajectory is approximated to be simple oscillatory motion z (t) ≈ z 0 cos (ωt) with angular frequency ω = 4 3νv s /H 5 1/2 and wavelength λ ≈ 2πv s /ω. Although µ = 0 in this study, a source dipole moment also leads to simple oscillations (Appendix C 2). This also theoretically explains the observations of persistent oscillations for neutral squirmers made by Zhu et al. [9], even though there are differences in the confining geometry. Next we add the dipolar term to the expansion where α = 3κ/H 3 , which is negative for pullers. The dipolar term also modifies the frequency ω due the wallinduced rotation Ω HI , but this effect is negligible if ν 81κ 2 / (48Hv s ), which is the case here. A pusher also obeys equation (3) but with α > 0 and exponentially growing amplitudes, which leads to a rapid breakdown of the near-centerline assumption. The sensitivity of oscillations to channel height is unmistakable in the exp H −3dependence of (3) reflecting the fact that the essential hydrodynamic moments are high order. Whereas higher order moments are required to predict the oscillation wavelength and damping factor quantitatively, the dipolar and quadrupolar moments can be fit from the dynamics using (3) with a error margin of ∼ 40% compared to the LBmeasured values. Hence, (1)-(3) allow for characterization of the swimmer's hydrodynamic properties based on experimental trajectories and can be readily transferred to the observations made by Zhu et al. [9]. Likewise, LB raspberry simulations can be extended to more complex 3D geometries such as square channels and round tubes, in which we observed helical motion (Appendix A 7).
In conclusion, we have studied the onset of oscillatory motion of swimmers in microchannels without externally applied flow and in an otherwise quiescent medium using both LB simulations and hydrodynamic theory. The pusher-type swimmers follow a sinusoidal trajectory with increasing amplitude, whereas pullers perform a damped oscillation towards the center of the channel. Our results and previous observations of such phenomena [9] can be explained by our theoretical model, which uses far-field hydrodynamics only. We conclude that the onset of oscillations can be described without taking into account near-wall lubrication effects as has been previously presumed [15] provided that a quadrupole moment (or source-dipole) is accounted for in addition to the pri-mary dipole moment. To fully characterize particle trajectories in relatively wide channels, many hydrodynamic moments are required, as high as the octupole in our case. However, the excellent match of our trajectory-fitted moments to those measured in bulk suggests that similar experimental measurements can be used to determine the hydrodynamic moment decomposition of microorganisms or artificial swimmers. Our work stresses the relevance of far-field hydrodynamics in confining geometries and thus opens the way for new studies that aim to exploit these insights in microfluidic environments. Future work will focus on the analysis of more complex force/counter-force swimmers to model the richness in shape and hydrodynamic moments of experimentally available swimmers.
Acknowledgements -AJTMM and TNS gratefully acknowledge funding from the ERC Advanced Grant In this section we present details of the simulation model that complement the description given in the main text. In addition, we provide further simulation results, which serve to underpin the generality of our findings.

Fluid Parameters and Swimmer Models
The 'raspberry swimmers' are based on the lattice-Boltzmann method implementation [20] and simulated using a graphics processing unit (GPU) based LB solver [24] that is attached to the MD software ESPResSo [22,25]. This GPU LB employs a D3Q19 lattice and a fluctuating multi-relaxation time (MRT) collision operator [26]. All of our simulations are performed in a quiescent (unthermalized) LB fluid. A threepoint interpolation stencil [27] is employed together with the LB viscous coupling of Ref. [28] to couple the raspberry particles to the fluid. We set the fluid density to ρ = 1.0m 0 σ −3 , the lattice spacing to 1.0σ, the time step to ∆t = 0.005τ (τ is the time and m 0 the mass unit), the (kinematic) viscosity to ν = 1.0σ 2 τ −1 , and the bare particle-fluid friction to ζ 0 = 25m 0 τ −1 . Fischer et al. [23] provide a detailed description of the dimensionless numbers that specify the fluid properties to which these choices correspond. We consider two types of self-propelled particles, a rod and cylinder as shown in Fig. 4. The rod consists of nine coupling points spaced 0.5σ apart over a line, with σ the LB grid spacing. The cylinder consists of 161 coupling points spread over 23 groups of hexagonal disks (seven particles with distance σ), stacked alternatingly with a separation of 0.5σ along the axis. The rod has an effective hydrodynamic length of 5.5 and a diameter of 1.7; the cylinder has an effective length of 12 and diameter of 3.2; and the sphere an effective radius of 3.1 [20]. Full details of these swimmers construction (mass, rotational inertia, etc.) are given in Ref. [20].
The raspberry bodies are made into swimmers by assigning a unit (direction) vectorp to their center that points along the symmetry axis, see Fig. 4a. Thisp comoves with the particle. We apply a force F in the direction ofp ( F = Fp) to the central molecular-dynamics bead, to which the rest of the coupling points are rigidly attached. This force causes the raspberry particle to move. We further apply a counter force − F to the fluid at a position lp, with l the separation length, to make the system force free. For positive values of l the swimmer is a puller and for negative values it is a pusher. We refer to Ref. [20] for the specific parameter choices. For convenience, we summarize the relevant quantities that these choices lead to in Table II, namely: the speed and hydrodynamic moments.
Using the size and speed of the swimmers, and kinematic viscosity of fluid, it is clear that Reynolds number of all our swimmers is less than 0.01. We use a quiescent fluid, therefore the Péclet number is ill-defined, as there is no translational (or rotational) diffusion. Both rod-and cylinder-type swimmers model 'cylindrical' self-propelled particles, but the level of description varies as well as the speed of the simulation. The rod-like model captures some of the hydrodynamic aspects of an extended object, namely the existence of a hydrodynamic quadrupole. The use of the low number of coupling points makes the   [20]. The table provides the shape, the velocity νs of the swimmer in units of (σ/τ ), the dipole strength κ (σ 3 /τ ), the quadrupole strength ν (σ 4 /τ ), the source-dipole strength µ (σ 4 /τ ), the source octupole o1 (σ 5 /τ ), and the force octupole o2 (σ 5 /τ ), respectively. The positive signs of κ correspond to pusher swimmers and the negative ones to pullers.
simulations fast compared to those for the cylinder swimmer. However, extended objects with a higher couplingpoint density, like our cylindrical swimmer, more accurately model a rod-like shape that is impenetrable to the fluid [29]. The use of a cylindrical swimmer thus serves to verify that the results obtained for the rod swimmer are not induced by low coupling-point density.

Simulation Setup
The above LB and raspberry coupling parameters result in faithful reproduction of theoretical results for passive particles in confining geometries, as was shown by De Graaf et al. [30]. Since we use a three-point coupling stencil deviations from the expected behavior of passive particles (solutions to the Stokes' equations) will occur within 2σ of the wall, rather than the 1σ found in Ref. [30]. Here, confinement is achieved by placing two no-slip (bounce-back) walls on either side of the simulation box in the z-direction. We pad the box using two lattice nodes of wall (zero velocity) on either side rather than one, because of the 3-point coupling of our swimmers to the fluid. The simulation domain is kept periodic in the other (xy) directions. This leads to a socalled 'slit-pore' geometry. We performed Poiseuille flow experiments to verify the height of the channel, the results of which we fit to the Hagen-Poiseuille expression. In all cases, the deviation between the imposed and fitted channel height is minimal (∼ 0.1σ).
In each of our simulations, the swimmer is initialized in the center of the box in the xy-direction, at a height z with respect to the center of the channel (z = 0). We ensure that the swimmer's directorp is in the xz-plane and impose its initial angle φ with the plate normalẑ, where φ ∈ [−π/2, π/2]. Typically, we use φ = 0 as the initial angle, which means that the swimmer is oriented parallel to the plane. The fluid in the channel is fully quiescent at time t = 0 and the particle starts with zero velocity. The LB parameters are chosen such that after the particle has moved only a fraction of σ the fluid flow field and terminal velocity of the swimmer is established, thereby minimizing the effect of inertia and momentumtransport retardation that physically do not play a role on the colloidal length scale at low Reynolds number.
To prevent the swimmers from penetrating the wall, we include a Weeks-Chandler-Anderson (WCA) interaction between the raspberry coupling points and the bounceback boundaries. The expression for the interaction is given by where r is the minimal distance between a coupling point and the wall and d is the 'diameter' of the particle. Every coupling point interacts with the wall via the WCA potential, leading to an overall wall-swimmer interaction that models that of a hard rod or cylinder with a hard wall. We typically use d = σ.

The Angular Evolution for Oscillating Swimmers
For completeness Fig. 5 shows the way the angle φ evolves along the trajectory of the swimmers given in Fig. 1 of the main text. The orientation of the rod changes along the trajectory. When the particle moves between the two walls, it comes close to making a 45 • angle with respect to the horizontal. This is further visualized in the supplemental movies described in the next section.

Description of the Supplemental Movies
To illustrate the movement of the swimmers, we have included two movies (not available in the arxiv version). These show the trajectory of a puller and pusher rod in a confining channel with height H = 13σ and lateral extent L = 70σ. The initial position is z = 1σ and φ = 0 and we used a WCA parameter of d = σ. The labeling of the movies is as follows: the type of the particle is given, followed by a list of quantities and values, with each set separated by a double underscore. The notation of the quantities is the one used throughout and each quantity and value are separated by a single underscore. We chose a slightly smaller lateral extent of the channel than used to produce Fig. 1 of the main text. The reason is that for the typical channel sizes studied in our work, the motion of the swimmer would be difficult to observe. However, we have verified that the limited size of the channel does not strongly effect the trajectory. Results are for swimmers that are initially oriented parallel to the walls and at z = 1σ; using the exact same data sets as used in Fig. 1 of the main text. (a) The results for pushers: rod for H = 13σ (red, solid), cylinder for H = 40σ (blue, dashed), and rod for H = 50σ (green, dots). (b) The results for pullers, otherwise the systems are the same.

LB Algorithm Limitations in the Near-wall Region
We scrutinize the presence of the artificial limit cycle for our pusher-type rod through a series of computational examinations. By our examinations we reach the following conclusions. Since the LB algorithm does not explicitly account for near-wall lubrication corrections [23], it fails to be accurate in the near wall regime and we are therefore unable to comment on the nature of any potential limit cycle. Additionally, the counter-force point can artificially penetrate the wall at the point of closest approach. These points indicate that, although limit cycles may exist in certain physical swimmers, the simulated trajectories cannot offer physically relevant predictions. We explain the way we arrived at these conclusions in detail below.
In our examination of the system, the lateral extent of the domain is varied between L = 5H and L = 35H to eliminate the effect of xy periodicity on our results: there is no discernible impact of L on the trajectories above L = 10H. We vary the viscosity and swimming speed to verify that retardation of the fluid momentum transport does not introduce these cycles; these changes only have a small effect. The value of the WCA interaction d is varied, as shown explicitly in Fig. 6. We find that for d > 1.5σ the limit cycle disappears and the rod's trajectory is reminiscent of the pusher cylinder's, see Fig. 1 in the main text. In both of these cases (inflated WCA rod and the unmodified cylinder) non-hydrodynamic contact with the WCA wall occurs and the self-propelled particles move along the plane of contact (sliding). Similar sliding dynamics have been observed in simulations that neglect HIs [31]. The pusher rod performs its persistent oscillatory trajectory even in the absence of the WCA potential. Fortuitously, it does not penetrate the wall, although penetration can be achieved in this case by starting with values of φ that are greater than ∼ 25 • when d = 0. This may seem to indicate that the limit cycle is a physical effect. However, this is not the case, as the rod comes very close to the wall, where LB does not faithfully reproduce hydrodynamics [30].
We therefore also considered the interaction of the counter-force point with the wall, see Fig. 7. We find that the counter-force interpolation (which takes place over a region of 3σ in diameter due to the three-point coupling) is partially inside of the wall at the closest approach, which impacts the reorientation of the rod. To check the effect of this, we switched to a two-point interpolation stencil. The limit cycle persists, but here too the interpolation region overlaps with the wall nodes, even though the overlap is substantially reduced. When the value of d increases beyond d = 1.5σ, the counter-force point is no longer interpolated inside the wall. Similarly, the cylinder's size prevents its counter-force point from being interpolated into the wall at closest approach. This indicates that the limit cycle observed for LB-raspberry swimmers is due to limitations in simulating the hydrodynamic interactions for close swimmer-wall separations, because of the spread-out counter-force scheme. While the near-wall hydrodynamics are not accurately captured by our algorithm, the far-field is. Therefore, in a system where there is a long-range (nonhydrodynamic) repulsion, our algorithm would produce the correct physics -provided that the range of the repulsion is sufficient to keep the LB coupling points far enough away from the wall. In the main text, we chose the WCA repulsion in such a way that the size of the 'hard core' matches the effective hydrodynamic size of the particle. Choosing the WCA range much larger, would remove this physical correspondence; therefore using an additional soft potential would be more appropriate to achieve wall repulsion. We are, however, unaware of any biological or artificial swimmers that are strongly repulsed from boundaries by long-ranged potentials and have therefore not considered this possibility further here.
In summary, the persistent oscillation (limit cycle) seen after long times for pusher-type swimmers must be attributed to a simulation artifact. Nevertheless, for the onset of the oscillation, which we are interested in the main text, there are no counter-force-overlap problems, as is clearly illustrated in Fig. 7.

Rod Swimmer Length
The effect of rod length on the trajectories of pushers is seen in Fig. 8. Since the effective hydrodynamic diameter of the rods is governed primarily by the coupling parameters when using only a single row of coupling points [29], varying the length has the effect of varying the aspect ratio of rod-shaped particles. We found only minor modifications of the trajectories, reflecting the minor changes in the hydrodynamic multipole expansion due to the change in aspect ratio. That is, the presence of a hydrodynamic quadrupole is the dominant effect in the formation of oscillatory trajectories; the strength of the quadrupole moment is only weakly perturbed by the changes in the length that we considered. The side length/diameter is H = 17σ. The swimmer is initially oriented parallel to the walls φ = 0 and located at x = 2σ and z = 1σ, with σ the MD unit of length and x = z = 0 the center of the tube. The trajectories shown here are 28H long in the y direction.

Helical Trajectories
The LB-raspberry swimmer model can be extended to simulate swimmers in other geometries. We find that our rods display helical trajectories, see Fig. 9. The helical trajectory observed for a tube with a square cross section is due to the swimmer starting off-center and away from one of the symmetry planes. Puller rods move consistently towards the center of the tube (not shown here) and also exhibit helical motion. The helical trajectory of the swimmer in the circular tube is due to a numerical artifact close to the boundary. The first part of the trajectory in the circular tube is purely oscillatory, as expected on the basis of symmetry and as we also observed for pullers. Only when an initial yaw angle is imposed does the rod perform a helical trajectory from the start of the simulation, similar to the observations of Ref. [9]. summed over. The two image systems due to point forces in the directionp are then p j B ± ij . From this pair of Stokeslet images, the image systems of the Stokes dipole, quadrupole, etc. can be constructed accordingly by taking successive derivatives as in equations (B3-B7) and the complete image system for the pair u ± is found.
These image velocity fields interact hydrodynamically with the swimmer. The wall-induced translational and rotational velocities of the force-free and torque-free swimmer are found by rearranging the Faxén relations evaluated at the swimmer position. Hence, we have where the derivatives act on the position x, E is the rateof-strain tensor, a is characteristic size of the swimmer, and G = γ 2 −1 γ 2 +1 is a function of the aspect ratio γ. Inserting the images of the swimmer-generated flow field (B1) into the Faxén relations (B10-B11) yields the wallinduced advection and rotation ( v HI and Ω HI in the equations of motion of the main text).

Appendix C: Swimmer dynamics model
Using the translational invariance along the x and y directions, we write the swimmer's orientation asp = (cos φ, 0, sin φ) without loss of generality, where φ = 0 corresponds to swimming parallel to the walls. Hence, the swimmer's equations of motion simplify to the two coupled equations,φ =φ(φ, z) andż =ż(φ, z). If we consider the simplified case of a point swimmer with aspect ratio γ = 1, and only use one image system on each side of the channel, these equations arė

Fitting Hydrodynamic Moments
An extension of these equations of motion (C1-C2), with a = 0 and G = 0, is used to match the dynamics of the model swimmers and LB swimmers (Table 1; main text). To achieve this, the time derivativesφ andż are extracted from the LB trajectories for a number of randomly chosen (φ, z) coordinate points, N = 500. Note that the first point in time is chosen to be after the first half oscillation such that the LB-raspberry swimmer has reached a constant swimming velocity and retardation effects are minimized. At each point, the LB values are compared to the values predicted by the model with a least squares method: Hence, the theory and LB simulations are matched by minimizing the S function with respect to the far-field multipole expansion parameters. Here, the swimming speed v s is fixed at the actual values (Table 1; main text). Likewise, the particle radius a is chosen to be fixed at the half-length of the LB rod or cylinder swimmer and the aspect ratio γ is set to its geometric value. Similarly, the parameters µ, o 1 are constrained to the LB-measured values, which is physically reasonable because these source doublets and quadrupolets are expected to be comparatively small, since our swimmers are constructed with-out fluid sources or sinks [20]. Finally, the multipole moments κ, ν, o 2 are allowed to vary, where a standard simulated annealing algorithm is used to find the least squares.

Analysis of swimmer oscillations
In order to analyze the micro-swimmer dynamics, we linearize the equations of motion (C1-C2) about the centerline of the micro-channel (z = 0), and about the orientation parallel to the walls (φ = 0). For simplicity, we consider only the dipolar and quadrupolar contributions to the multipole expansion and set the octupolar and higher-order contributions to zero. The dynamics can then be captured by the matrix equation First, we consider the motion in the absence of a dipole moment (κ = 0), but with quadrupole moment ν and source doublet moment µ. Then, the eigenvalues λ e of the matrix are λ e = ± 4 3(ν + 2µ) (14ν + 16µ − v s H 3 )