Squirmer hydrodynamics near a periodic surface topography

The behaviour of microscopic swimmers has previously been explored near large-scale confining geometries and in the presence of very small-scale surface roughness. Here, we consider an intermediate case of how a simple microswimmer, the tangential spherical squirmer, behaves adjacent to singly and doubly periodic sinusoidal surface topographies that spatially oscillate with an amplitude that is an order of magnitude less than the swimmer size and wavelengths that are also within an order of magnitude of this scale. The nearest neighbour regularised Stokeslet method is used for numerical explorations after validating its accuracy for a spherical tangential squirmer that swims stably near a flat surface. The same squirmer is then introduced to different surface topographies. The key governing factor in the resulting swimming behaviour is the size of the squirmer relative to the surface topography wavelength. For instance, directional guidance is not observed when the squirmer is much larger, or much smaller, than the surface topography wavelength. In contrast, once the squirmer size is on the scale of the topography wavelength, limited guidance is possible, often with local capture in the topography troughs. However, complex dynamics can also emerge, especially when the initial configuration is not close to alignment along topography troughs or above topography crests. In contrast to sensitivity in alignment and topography wavelength, reductions in the amplitude of the surface topography or variations in the shape of the periodic surface topography do not have extensive impacts on the squirmer behaviour. Our findings more generally highlight that the numerical framework provides an essential basis to elucidate how swimmers may be guided by surface topography.


Introduction
Ever since the studies of Robert Hooke and Antonie van Leeuwenhoek, it has been known that a drop of pond water contains countless microbes, many of which are motile, and indeed some can be lethal, such as Naegleria fowleri, the causative agent of amoebic meningitis. Even harmless motile microbes have attracted the attention of scientists for centuries, though more recently developments in nano-and micro-technology have also enabled fabrication of self-propelling artificial micro-robots and manipulation of their dynamics using microfluidic devices (Kherzi and Pumera, 2016). In laboratory experiments and observations, with both synthetic and biological swimmers, of the range of known control mechanisms, by far the most common is the influence of confining boundaries.
These curved boundaries and obstacles are typically larger than, or comparable to, the swimmer. If the structure on the wall boundary is smaller than the swimmer length scale, it may be considered a rough boundary instead of a completely flat surface. The impact of surface roughness has previously been considered via an asymptotically small amplitude of the surface topography in the presence of a spherical particle and a spherical microswimmer (Rad and Najafi, 2010;Assoudi et al., 2018), the so-called squirmer (Shaik and Ardekani, 2017).
The squirmer is a model microswimmer first proposed by Lighthill (1952) as a slightly deforming sphere and later corrected and used by Blake (1971) as a model of ciliate swimmers. This simple model is currently understood to provide qualitative predictions for a spherical biological microswimmer (Pedley, 2016;Pedley et al., 2016). In particular, a simplified version of the model, in which a rigid sphere can self-propel due to a given surface velocity slip profile, is known as the spherical tangential squirmer. This has been widely used as a simple mathematical model having a finite volume for studies on hydrodynamical aspects of microswimming such as nutrition uptake (Magar et al., 2003), cell-cell interactions (Ishikawa et al., 2006;Drescher et al., 2009), Janus particle motility (Spagnolie and Lauga, 2012;Ishimoto and Gaffney, 2013), collective dynamics (Evans et al., 2011;Zöttl and Stark, 2012;Delfau et al., 2016;Oyama et al., 2016), swimming in a non-Newtonian medium (Lauga, 2009;Zhu et al., 2012;Nganguia and Pak, 2018), and swimming near a wall (Llopis and Pagonabarraga, 2010;Spagnolie and Lauga, 2012;Ishimoto and Gaffney, 2013). The squirmer has also been studied in the context of confinement and obstacles such as the interior of a tube (Zhu et al., 2013), the presence of lattice-like multiple obstacles (Chamolly et al., 2017), or a curved and structured wall (Das and Cacciuto, 2019).
Investigations into the effects of small surface topography on microswimmers are, however, limited to the asymptotic analysis of rough surfaces or boundary features (Kurzthaler and Stone, 2021) such as curvatures with length scales that are much larger than those of the swimmer. The current study, therefore, aimed to bridge the gap between an asymptotically small amplitude of the surface roughness and large length scale curved boundaries. For periodic structures at this mesoscale, in particular, there is the prospect that the microswimmer may become oriented and guided by the surface, and we will numerically investigate the dynamics of a spherical tangential squirmer under these conditions. Such investigations are particularly motivated by recent studies of a colloidal microswimmer near a small surface topography (Simmchen et al., 2016), highlighting that a structured surface topography may be fabricated in a microfluidic device with the potential for utilisation in guiding microswimmers.
The very near-wall dynamics, at a separation of around 50 nm and less (Klein et al., 2003), typically depends on both hydrodynamic and steric interactions (Klein et al., 2003;Kantsler et al., 2013;Bianchi et al., 2017), and a short-range repulsive potential force is often utilised by modelling studies to ensure that simulated cells do not penetrate the walls (Spagnolie and Lauga, 2012). However, even a small difference in the repulsion function can alter swimmer behaviour (Lintuvuori et al., 2016;Ishimoto, 2017). Thus, in this initial study, we only focus on the hydrodynamic interactions and do not consider any additional steric interactions, contact mechanics, and charge effects.
This is motivated not only by the utility of the relative simplicity in this initial study but also for understanding the impact of hydrodynamic surface interactions where, despite these non-hydrodynamic forces, swimmer deposition is undesirable and thus of minimal interest, compared to topography guidance dynamics when deposition does not occur. It also entails that the results and conclusions of this study are not contingent on the details of contact mechanics and steric forces, which vary with surfaces and solutes (Klein et al., 2003). Thus, the numerical simulations are stopped just before the squirmer dynamics is influenced by the short-range dynamics on a very close surface approach. This short-range detail may be accommodated in later work together with many further refinements, such as incorporating more faithful representations of flagellated and ciliated swimmers (Ohmura et al., 2018;Ohmura et al., 2021).
The structure of this paper is as follows: Section 2 introduces the squirmer model and three different surface topographies. Section 3 discusses the numerical methods and their verifications. Section 4 presents the simulation results for the different surface topographies, followed by a discussion of the implications, in particular for microswimmer guidance via surface topography, in Section 5.

The squirmer
We consider the non-dimensional Stokes equations of the low Reynolds number flow, from which it follows for a velocity field u and a pressure field p that ∇p Δu and ∇ · u 0. (1) We impose the no-slip boundary condition on the swimmer and the wall, together with the force and torque balance equations for a swimmer with negligible inertia.
Frontiers in Cell and Developmental Biology frontiersin.org 02 We first introduce the spherical tangential squirmer model. The no-slip boundary is imposed on a spherical swimmer of dimensionless radius a = 1, possessing an axisymmetric and tangential surface velocity (Ishikawa et al., 2006). The sphere centre is located at X = (x, y, z), and n denotes the unit vector of its orientation, as shown schematically in Figure 1, where coordinates, variables for the position of the squirmer centre, and the angles θ, Θ are defined using a diagram. Here, θ is the angle made by the swimming direction and the xy plane. In particular, Θ ∈ [0, π] denotes the polar angle relative to n, and we impose an axisymmetric tangential velocity slip u s on the squirmer, given by where V n is a function derived from the Legendre polynomial P n (x) via The swimming velocity in free space is dictated by the first mode, with U = (2/3)B 1 n (Lighthill, 1952;Blake, 1971). We fix B 1 = 3/2 so that the squirmer swimming speed is set to be unity in free space, and we neglect the higher modes by setting B n = 0 for n ≥ 3 so that the swimmer is subsequently fully characterised by the squirmer parameter β = B 2 /B 1 (Ishikawa et al., 2006). In particular, and following convention, the swimmer is denoted as a pusher when β < 0, a puller when β > 0, and a neutral swimmer when β = 0 (Evans et al., 2011). The second mode, associated with the parameter B 2 , corresponds to the flow induced by the Stokes dipole. In particular, a cell with a trailing flagellum, such as an E. coli bacterium or a sperm cell, behaves as a pusher; cells with leading flagella, such as Chlamydomonas and Leishmania (Walker et al., 2019), are modelled as pullers; and cells possessing fore-aft symmetry, such as ciliates, behave as neutral swimmers (Evans et al., 2011).
Here, we focus on spherical tangential squirmers that swim stably near a flat surface. Thus, we consider puller squirmers with β ≳ 5, which are known to exhibit stable swimming near a flat surface (Ishimoto and Gaffney, 2013;Uspal et al., 2015). In particular, we examine their dynamics close to surfaces with structured topographies. The first of these topographies is a onedimensional sinusoid defined by where A is the amplitude and k = 2π/λ is the wavenumber, with λ denoting the wavelength ( Figure 2A). The second topography is given by the doubly periodic sinusoid.
as depicted in Figure 2B, and the third is given by This final topography is shown in Figure 2C and presents doubly periodic peaks with highest and lowest heights of +A and −A, respectively, as in the previous two topographies. However, notably, the inter-peak wavelength is now halved, and the parameter λ = 2π/k no longer represents the wavelength since the sinusoidal functions are squared in Eq. 6. Throughout this study of the doubly periodic topographies, we focus on cases that are symmetrical in switching the x and y-directions. We consider both the surface of the squirmer and the wall surface topography denoted by S and W, respectively, with the boundary conditions of the Stokes equations given by no-slip conditions on both boundaries. The surface velocity of the squirmer, denoted v(x), can be described by combining the squirmer linear velocity U and angular velocity Ω, together with its tangential surface velocity, u s , of size given by Eq. 2 and in the axisymmetric tangential direction, as depicted in Figure 1. Hence, the no-slip condition entails the fluid velocity on the surface of the swimmer is given by In contrast, the wall surface topography is assumed to be stationary, and we thus have 3 Numerical methods

Nearest-neighbour regularised Stokeslet method
The dynamics of the squirmer has been computed using the nearest-neighbour regularised Stokeslet method (nnRSM)

FIGURE 1
Schematic diagram of a spherical tangential squirmer, of radius a = 1, near a no-slip wall. Here, X denotes the centre of the spherical swimmer relative to a laboratory reference frame, with Cartesian coordinates (x, y, z) and h(x, y) is the height of the surface above its average midplane (dashed) at z = 0. Furthermore, z is overloaded and also represents the height of the swimmer centre above the midplane, with an analogous overloading of x, y. Whether x, y, z refer to coordinates or the overloaded definition X = (x, y, z) for the location of the squirmer centre will be clear from context. The unit vector n gives the orientation of the swimmer, which makes an angle θ relative to the mid-plane of the surface topography, and Θ is the local polar angle of a point on the swimmer' surface relative to n. The swimmer's motility is driven by an axisymmetric tangential velocity of its surface, of size u s , and in the direction of increasing Θ, as detailed in the main text.
Frontiers in Cell and Developmental Biology frontiersin.org 03 (Gallagher and Smith, 2018;Smith, 2018), and the numerical simulations have been performed based on the MATLAB code accompanied by Gallagher and Smith (2018), as we now summarise. The Stokes flow boundary integral equations for the single-layer formulation are given by (Pozrikidis, 1992) Here, f i (y) denotes the components of the surface traction and the integral kernel S ij is the Stokeslet, which exhibits an integrable singularity as y → x, with the surface integral well-defined. For numerical tractability, Cortez introduced a regularised Stokeslet (Cortez, 2001), which is the exact divergence-free solution to the Stokes flow equations with a spatially smoothed point force, and then approximated the boundary integral (Cortez et al., 2005) via where S ϵ ij is the regularised Stokeslet and ϵ is the regularisation parameter. As the error in this approximation is O(ϵ), we recover the singular boundary integral once we take the limit of ϵ → 0. Following Cortez et al. (2005), we consider a regularised Stokeslet of the form, where r = x − y, r = |r|, and δ ij represent Kronecker's delta. The noslip boundary conditions are simply given by enforcing u(x) = v(x) from Eqs. 7, 8 for boundary points in the integral Eq. 10. Since the squirmer is swimming freely, we also have the inertialess force and torque balance equations.
We then discretise the surface integrals (10), (12), by introducing the quadrature node positions x[n] and the associated weights A[n] for the discretised surface point n (n = 1, . . ., N), where N is the total number of surface points. The aforementioned surface integrals contain the product of 'f dS' and we discretise the integral (Gallagher and Smith, 2018;Smith, 2018) via where the symbol, •, on the right-hand side, represents an arbitrary Continuing with the framework of Smith (2018), we introduce a second surface discretisation, x[q], (q = 1, . . ., Q) which corresponds to a more refined discretisation than used for the surface traction, with N ≪ Q, as illustrated in Figure 3. The regularisation error is O(ϵ), and this motivates taking relatively small values of ϵ in computations. The two discretisations enable an efficient numerical solution as the kernel, S ϵ ij , which can vary rapidly, and thus requires a finer discretisation, which would be inefficient if used for the surface traction, f. In particular, the size of the dense linear system depends only on N, not Q. Thus, the cost of assembling the system is defined by O(NQ), whereas the cost of the direct solver is defined by O(N 3 ). Moreover, if the force and quadrature discretisations do not overlap, the quadrature error no longer diverges as ϵ → 0, and hence a less refined force discretisation in this framework is in general more accurate than if the two discretisations coincide (Gallagher et al., 2019).
The nearest-neighbour matrix, ][q, n], is then defined separately for a swimmer and a wall as where arg min is the argument at which the minimum is attained over the setn ∈ {1, . . . , N}, and we use this matrix to interpolate the discretisation via Combining Eqs 13, 15 and noting the total number of the points for the finer discretisation Q, we have

FIGURE 2
Illustrations of the surface topography functions h(x, y) considered in this study. (A) Singly periodic sinusoidal wave topography, (B) doubly-periodic sinusoidal wave, and (C) doubly-periodic peaks, respectively, associated with Eqs 4-6. The maximum and minimum heights, h(x, y), are +A and −A for all topographies.
In particular, to solve the squirmer trajectory, we first need to determine its velocity, U, and angular velocity, Ω, which can then be integrated over time to determine the squirmer location and orientation. First, it should be noted that at a fixed point in time, the squirmer location and orientation are obtained from previous integration or from the initial conditions at the start of the simulation. Then, the discretisations of Eqs 10, 12, with u = v eliminated in terms of U, Ω, and the known u s via Eqs. 2, 3, 7, 8, give 3N + 6 constraints for the 3N + 6 scalar unknowns associated with the unknown surface tractions at the N discretisation points and the unknowns U, Ω. The resulting linear system is readily solved, provided that collocation points are unique, and we have a nonsingular dense linear system that can be solved directly via standard methods.
As is the case for both singular and regularised versions of the boundary integral representation for flow around a constant volume body, the integral equation admits a gauge freedom f → f + αm, where α is any constant and m is the surface normal pointing into the fluid (this can be observed by applying incompressibility and the divergence theorem to deduce that S S ij m i dS = 0). In the absence of boundary conditions for traction, this freedom results in the pressure being determined only up to an additive constant in the exact problem. Discretisation results in an invertible matrix, and hence a unique (approximate) solution, because the discretised integral is no longer evaluated precisely to 0; moreover, the nonuniqueness of the continuum solution for the pressure is not dynamically important as it does not affect either the total force or moment on the swimmer.

Swimming in a free space
We first examine the numerical accuracy of the swimming velocity calculation for the squirmer in free space. The squirmer parameter is set to β = 0, and the exact swimming speed is |U| = 1, as detailed in the previous section. We have fixed the regularisation parameter ϵ = 0.001 and examined the impact of changing the discretisation refinement. In particular, with the total number of the points that form the squirmer surface given by N 6n 2 s and Q 6N 2 s for each discretisation (Gallagher and Smith, 2018;Smith, 2018), changes in both n s and N s have been examined. The results of Table 1 establish numerical parameters sufficient to obtain our desired relative error of around 1%. Finally, we also note that changes in β have not been observed to alter the swimming speed, as expected.

Swimming near a wall
Hereafter, we set the squirmer discretisation parameters to be (n s , N s ) = (4, 18) and consider the squirmer near a no-slip wall. As previously studied by the boundary element method (Ishimoto and Gaffney, 2013), a strong puller tends to stably swim near a flat wall. We, therefore, choose β = 7 and set the initial location of the squirmer centre to be (0,0,1.15), with the initial orientation given by θ = −0.17π, which is effectively the initial angle of attack relative to the mid-plane of the surface topography, as can be seen from Figure 1. Illustration of the points representing the squirmer and surface topography, with the finer discretisation (A) used for the kernel and the coarser discretisation (B), used for the surface traction. Discretisation points on the spherical squirmer surface are indicated by red dots. In contrast, blue discretisation points are shown for the doubly periodic surface topography of Figure 2B, with the size of the discretised surface given by L = 8 and the surface topography given by Eq. 5, as plotted in Figure 2B, with an amplitude of A = 0.1 and the wavelength given by λ = 2. Note that we here display the wall meshes after rescaling by the method described at the end of Section 3. We first use the regularised Blakelet (Ainley et al., 2008) in the nearest-neighbour regularised Stokeslet method and compare the simulation result with the stable distance obtained via the boundary element method using the singular Blakelet (Ishimoto and Gaffney, 2013), with the latter providing the stable separation distance z* ≈ 1.1578. As shown in Figure 4, these predictions of these two algorithms are in reasonable agreement.
We then consider a wall that is captured by an explicit discretisation of its surface rather than by use of the Blakelet, as implemented in the sperm simulation of Gallagher and Smith (2018). The wall is given by the x-y plane and represented by the square with a length of L, with its centre at the location of the projection of the squirmer centre, X, onto the plane z = 0. Each side contains n w and N w points with equal separations for surface traction and kernel discretisations, respectively (Gallagher and Smith, 2018;Smith, 2018). Hence, the number of points on the surface S ∪ W are given by N 6n 2 s + n 2 w and Q 6N 2 s + N 2 w for the surface traction and kernel discretisations, respectively. An example of a swimming trajectory is plotted in Figure 4.
We then rescale the wall points to resolve the squirmer-boundary hydrodynamic interaction more efficiently by using the function The equally discretised square of unit length S x, y ∈ −1 2, 1 2 × −1 2, 1 2 is mapped by this function, via and then dilated by the scale of L. The square obtained by this scheme more precisely represents the hydrodynamical interactions between the squirmer and the wall, as seen from the results labelled "rescaled" in Figure 4, which use this mapping. These trajectories in particular are sufficiently accurate for our purposes and very close to the prediction of the boundary element method (BEM) of the stable swimming height above the surface, which is exact to within discretisation error.

Results
In this section, we discuss the swimming trajectories of the squirmer near a surface with a structured periodic topography, as defined in Eqs. 4-6 and depicted in Figure 2. For all simulations presented, the initial height of the squirmer was fixed at z = 1.2, with the initial angle of attack given by θ = −0.17π. Furthermore, initial squirmer centre location coordinates of x = 0, y = 0, are set together with squirmer parameters of B 1 = 3/2, β = 7, and a surface topography amplitude of A = 0.1, unless explicitly stated otherwise. Although the dynamics can change depending on the parameters and initial settings, we fix these variables to consider the stable behaviour and its modulation by surface topography, focussing on the impact of the initial orientation and topographic patterns. The surface topography wavelength and the initial orientation of the squirmer in the x-y plane, namely, φ in Figure 5, are thus varied extensively among the simulations and either stated or, in the case of φ, can otherwise be immediately inferred from the initial tangent angles of the trajectories in the presented plots.

FIGURE 4
Swimming trajectories of the squirmer with a regularisation of ϵ = 0.001 and different discretisation parameters. Also plotted is the height of the stable fixed point obtained by the boundary element method using the singular Blakelet (Ishimoto and Gaffney, 2013), as labelled "BEM stable distance," and a trajectory using a regularised Blakelet with the nearest-neighbour regularised Stokeslet method, as labelled by "Blakelet."

FIGURE 5
Bird's eye view of the squirmer above the doubly periodic surface topography of Eq. 5, as depicted in Figure 2B, with λ = 4. The angle φ is defined to be the angle of the x-y projection of the squirmer orientation vector n from the y-axis, as illustrated, and thus gives the direction of the squirmer in the x-y plane.
Frontiers in Cell and Developmental Biology frontiersin.org

Singly periodic sinusoidal topography
We start with the single-wave sinusoid topography given in Eq. 4 and Figure 2A. The initial location of the squirmer and the initial angle of attack are fixed at the simulation default initial values of X = (0, 0, 1.2) and θ = −0.17π, as stated previously, while the initial orientation of the squirmer in the x-y plane has been initially considered in detail with φ = 0.5π and then subsequently varied.
Thus, we first consider a squirmer swimming parallel to the wave vector of the sinusoidal topography or equivalently along the x-axis. Fixing the initial orientation relative to the x-y plane via φ = 0.5π, swimming trajectories in the x-z plane are plotted in Figure 6A for different surface topography wavelengths λ = 1, 2, 4, 8. Corresponding trajectories in the θ-z phase plane are shown in Figure 6B. When the wavenumber λ is smaller (λ = 1, 2), the squirmer attains stable oscillatory swimming, but the oscillation in the z-axis is smaller than the surface topography amplitude, A = 0.1, highlighting that the topography only perturbs the stable position associated with swimming near a flat wall. This may be observed in Figure 6C, where the z-dynamics for the last part of the oscillating motion obtained in Figure 6A are shown relative to a horizontally rescaled and shifted axis x/λ together with the surface topography function h. However, as the wavelength is increased to λ = 4, λ = 8, the oscillatory motion then transitions to an amplitude that is larger than that of the surface topography, as can be seen in Figure 6C. In addition, one can observe that, with λ = 4, the wavelength of the z-component oscillations in the trajectory need not match that of the underlying surface topography, though in contrast, these two wavelengths do match for the trajectories with λ = 2 and λ = 8. Although the problem converges to the locally flat wall case in the large wavelength limit, the locally flat approximation does not hold for the range of λ we have examined, as may be inferred from the oscillatory motion in Figure 6C.
We then consider the squirmer dynamics with different initial values of φ, and thus different initial orientations relative to the x-y plane, while once again varying the wavelength of the singly periodic topography. Figure 7 shows the predicted trajectories and the orientations for the case with the surface topography amplitude A = 0.1 and wavelength λ = 1, 2, 4, 8, while considering various values of φ, from 0 to 0.5π. From the figure, one can observe that the squirmer tends to swim either parallel or perpendicular to the direction of the well, which is aligned along the y-axis, though the orientation angle φ need not always necessarily match the direction of motion. For instance, some trajectories in Figure 7 follow the topographical crest by swimming in the direction of the positive y-axis, with the orientation angle, φ, remaining very close to the initial value rather than aligning with the y-axis, even after a long time. Also, some further trajectories are attracted towards the negative x-axis, without the orientation angle φ evolving to reflect this change in swimming direction. Hence, overall drifting, that is, a misalignment between the swimming direction and swimmer orientation, can be induced by the squirmer-topography hydrodynamic interaction.
Notably, the swimming dynamics associated with an initial orientation angle of φ ≈ 0.05 with λ = 1, or φ ≈ 0.15 with λ = 2, is unstable, and the trajectories evolve towards the stable orientations of φ = 0, as seen in Figures 7A, B. Here, the stable swimming along the y-axis is accompanied by hydrodynamic Frontiers in Cell and Developmental Biology frontiersin.org 07 capturing, with the squirmer moving along a trough of the surface topography. Furthermore, in both cases, the z-dynamics attains a stable oscillatory motion, as may be observed in Figures 6A, B. In both cases, the qualitative aspects of these features are unaltered when the amplitude A is decreased, though the timescale for the reorientation along one of the axes becomes longer.
Analogously, an increase in the topography wavelength to the intermediate value of λ = 4 entails that swimming oblique to the troughs and peaks of topography can be observed, as seen in Figure 7C. However, at this wavelength, the squirmer concomitantly undergoes extensive oscillations in the z-direction. Furthermore, the squirmer enters the near vicinity of the surface ( Figure 7C) once it is no longer oriented approximately along the yaxis. This requires a consideration of surface mechanics to proceed. More precisely, in practical applications, surface mechanics would become relevant and would need to be added to the model. This is outside the detailed scope of the study and, hence, we stop the numerical simulation, thus also avoiding the numerically unreasonable spatial resolutions required for the associated finescale hydrodynamics.
We further increase the wavelength of the sinusoidal topography to λ = 8, with trajectories presented in Figure 7D, which on projection to the x-y plane, are essentially straight, regardless of the initial orientation angle φ. Hence, at larger wavelengths, the squirmer swims with a direction that is unaffected by the surface. Furthermore, the z-dynamics of the squirmer trajectory become more oscillatory as the topography wavelength increases, as observed previously in Figure 6, unless the squirmer is captured in the trough along the y-axis, in which case the z-dynamics converges to a stable position.
More generally, all of these observations highlight that even with a surface topography amplitude of A = 0.1, which barely visible, as highlighted by Figure 3, the squirmer's behaviour is affected by the structured surface topography in a complex manner. In particular, the resulting trajectories are contingent on the details of the topography parameters and the squirmer orientation, especially once the topography wavelength is comparable to the squirmer size.

Doubly periodic sinusoidal topography
We now consider the squirmer dynamics near a surface with the doubly periodic wave topography, given by Eq. 5 and illustrated in Figure 2B. In the current setting, the topography breaks the Frontiers in Cell and Developmental Biology frontiersin.org 08 translational symmetry in the y-direction, and hence the trajectories now also depend on the value of the initial y position of the squirmer centre. We first consider the squirmer starting with an orientation φ = 0.5π and an initial centre location coordinate of y = λ/4, with the other initial location coordinates and the initial angle of attack at the default values of x = 0, z = 1.2, and θ = −0.17π. Then, the squirmer moves along the −x-axis with similar dynamics in the z-direction to that displayed in Figure 6. In particular, the dynamics in the zdirection are only slightly perturbed when λ = 1, 2 but display larger oscillations when λ = 4, 8.
We then consider variations in the initial squirmer orientation angle φ that, as previously, entails the trajectories are no longer constrained to two spatial dimensions. The surface topography amplitude remains fixed at A = 0.1, and the initial height, location, and attack angle are at default values, while we consider variation in the orientation angle φ and the topography wavelength, which here is still given by λ. When λ = 1 ( Figure 8A), the trajectories are not affected by the surface topography, as the trajectories are straight when projected on to the x-y plane, and the angle φ is constant in time. However, in the simulation with the wavelength λ = 2 ( Figure 8B), some trajectories with φ ≈ π/4 are hydrodynamically captured near the bottom of the doubly periodic sinusoidal valley, whereas swimming outside of this region of initial orientation angles is not affected by the surface topography, as in the case of λ = 1.
These features can also be observed when we increase the wavelength to λ = 4 ( Figure 8C). In contrast, for the larger wavelength of λ = 8, there is no evidence for an attracting trough of squirmer dynamics near the bottom of the topographic valley ( Figure 8D). Together, these results highlight that the hydrodynamic attraction towards, and subsequently along, topographic valleys is not only limited but also only possible when the scale of the swimmer's diameter is comparable with the length scale of the surface topography.
We then move to consider the final surface topography of doubly periodic peaks as introduced by Eq. 6 and displayed in Figure 2C. Trajectories with straight line projections onto the x-y plane can be observed when the squirmer's initial location and orientation angle align along topography troughs or across topography crests. For example, given the default initial location of X = (0, 0, 1.2) and an initial orientation angle φ = 0, the squirmer swims with φ = 0 throughout time. In addition, these trajectories also exhibit nearly constant z-dynamics, though small z-oscillations are observed with amplitude ≲ 0.05 due to the topography, and this oscillation is further diminished as the topography wavelength increases. Furthermore, with the default initial location and an initial value of φ = 0.25π, or with initial values of X = (0, λ/4, 1.2) and φ = 0.5π, straight line x-y projected trajectories are also observed. Furthermore, in the z-direction, the squirmer behaviour changes from small oscillatory perturbations to larger amplitude zoscillations as the wavelength increases, in direct analogy to the examples considered in detail with the previous topographies.
For more general initial configurations of the squirmer, we have considered the three-dimensional behaviour of the resulting trajectories, with the simulation results plotted in Figure 9 for λ = 2, 4, 8, 16, observing the wavelength for this topography is λ/2, not λ. In the previous doubly periodic surface topography, swimming with the orientation angle φ = π/4 allowed the cell to move along a surface topography trough, noting that the prospect of drifting observed in the singly periodic topography of Figure 2A can be ruled out. In contrast, for the current case, an orientation of φ = π/ 4 moves across the surface topography peaks.
Again, the simulation results show that the squirmer is attracted by orientation angles corresponding to troughs in the surface

FIGURE 8
Dynamics of the squirmer near a surface with the doubly periodic sinusoidal wave topography of Eq. 5, as depicted in Figure 2B. The surface topography amplitude is given by A = 0.1, and the wavelength is (A) λ = 1, (B) λ = 2, (C) λ = 4, (D) λ = 8. (Top panels) the projections of the squirmer trajectories onto the x-y plane with different initial orientation angles, φ, as defined via Figure 5. These initial angles may be inferred from the initial tangents of the plotted projected trajectories. (Bottom panels) the time evolution of the orientation angle φ. Different colours index different initial conditions.
Frontiers in Cell and Developmental Biology frontiersin.org 09 topography, which here are φ = 0, π/2, for example. Furthermore, such topographic attraction is realised when the squirmer diameter is close to the characteristic length scale of the surface modulation, as seen in Figures 9A-C. However, again for the parameter regimes of Figures 9B, C, the swimmer can approach very close to the surface, and the trajectory simulations are halted as detailed surface dynamics are outside the scope of the study. Finally, we note that once the surface topography wavelength is increased further, with λ = 16 as presented in Figure 9D, all projected trajectories are straight lines, and the orientation angle φ is constant in time. Hence, at these larger wavelengths, the surface-swimmer interactions no longer influence the guidance of the squirmer.

Discussion
In this paper, we have numerically investigated the hydrodynamics of a puller spherical tangential squirmer near a surface with a singly or doubly periodic structured topography. In particular, the amplitude of the surface topography was an order of magnitude less than the squirmer size and a wavelength on the scale of the squirmer size. The simulated squirmer is known to be attracted to a stable separation from a flat wall, and a mesh-free regularised Stokeslet boundary element numerical scheme was demonstrated to accurately capture the dynamics induced by the subtle hydrodynamic interactions between a spherical tangential squirmer and a flat wall.
When the wavelength of the sinusoidal surface topography is smaller than the squirmer size, the perpendicular dynamics of the swimmer trajectory is a small amplitude oscillatory perturbation from the constant stable swimming height associated with a flat boundary. However, as the wavelength of the surface topography is increased, the squirmer acquires larger vertical, z-direction, oscillations with a wavelength that matches that of the topography at very large surface topography wavelengths, but not always at intermediate values ( Figure 6C).
Furthermore, the squirmer movement in the horizontal, x-y, plane has been observed to be highly dependent on the detailed geometrical properties of the surface topography. We first considered a singly periodic sinusoidal surface topography. When the wavelength of the surface topography (Eq. 4) is not significantly larger than the squirmer diameter (λ < 4), the horizontal squirmer motion reorientates towards one of two stable directions, i.e., parallel and perpendicular to the wavevector of the surface topography, as seen in Figure 7. Furthermore, drifting can sometimes be observed, whereby the direction of motion differs slightly from the orientation angle, φ, as noted in Figure 7, for example.
At intermediate wavelengths (λ = 4) the squirmer can approach the surface. The detailed subsequent behaviour would be contingent on the near-surface physics, the detailed study of which is beyond the scope of this study. Once the surface topography wavelength is further increased, with λ = 8 sufficient so that the wavelength is four times that of the squirmer diameter, we observed that the horizontal motion is that of straight lines and surface induced guidance of the squirmer in the horizontal plane is lost.
For the doubly periodic surface topographies, the squirmer had the tendency to be locally guided to swim along surface topography

FIGURE 9
Dynamics of the squirmer near a surface with the doubly periodic sinusoidal wave topography of Eq. 6, as depicted in Figure 2C. The surface topography amplitude is given by A = 0.1, and the wavelength is unchanged from previous plots, but is no longer given by λ = 2π/k where k is the wavenumber of the sinusoidal functions in Eq. 6, since these functions are squared. Hence, to preserve wavelength at 1, 2, 4, 8 length units in the respective columns, we take (A) λ = 2, (B) λ = 4, (C) λ = 8, and (D) λ = 16 on moving from left to right across the figure. (Top panels) the projections of the squirmer trajectories onto the x-y plane with different initial orientation angles, φ, as defined via Figure 5. These initial angles may be inferred from the initial tangents of the plotted projected trajectories. (Bottom panels) the time evolution of the orientation angle φ. Different colours index different initial conditions.
Frontiers in Cell and Developmental Biology frontiersin.org troughs rather than over crests, though only when the squirmer diameter was comparable to the surface topography wavelength. However, such tendencies were weak and, at intermediate surface topography length scales, often accompanied by very close approaches to the surface where steric interactions would be important. These qualitative features of the squirmer trajectories have also been observed when the squirmer parameter, β = B 2 /B 1, from Eq. 2, was varied within the range that induced stable swimming in the vicinity of a flat wall. In addition, changing the sinusoidal surface topography to a sigmoidal topography to represent fabricated surface wells in microfluidic devices did not induce a significant change in the qualitative features of the swimming trajectories. In turn, this evidences that squirmer swimming behaviours are influenced mainly by the length scale of the surface topography. In addition, we have also observed that the swimmer behaviour can be complex, especially when the swimmer is not aligned along the surface topography troughs or above the surface topography crests, though a limited local guidance to swim along the surface topography troughs has also been common.
There are considerable numbers of studies focussing on microswimmer dynamics near non-trivial geometrical structures such as curved obstacles (Nishiguchi et al., 2018;Das and Cacciuto, 2019), bumps (Simmchen et al., 2016;Yang et al., 2019), and mazelike micro-devices (Denissenko et al., 2012;Tung et al., 2014), but the length scale of the surface topography in the current study features much finer surface structures. We also note that the transitional vertical behaviours in the z-direction, from perturbations of the stable swimming height for a flat wall to topography-following motion at very large surface topography wavelength, necessitate a consideration of the finite-size amplitude of the surface topography. In particular, such behaviours highlight that the dynamics examined in this study requires larger-scale physics beyond the effective boundary conditions (Sarkar and Prosperetti, 1996;Kunert et al., 2010;Luchini, 2013) based on a very small amplitude surface roughness.
In the context of representing biological microswimmers such as spermatozoa and bacteria, which are pusher swimmers rather than pullers, hydrodynamic stable swimming occurs for prolate pusher tangential squirmers, but not spherical squirmers (Ishimoto and Gaffney, 2013). Moreover, the hydrodynamic interactions strongly depend on the swimmer morphology and beating pattern (Smith et al., 2009;Shum et al., 2010;Ishimoto and Gaffney, 2014;Walker et al., 2019). In turn, this highlights that detailed numerical studies are required to explore the surface dynamics for both prolate squirmers and more realistic microbiological swimmers near non-trivial surface topographies, including the prospect of a ciliated epithelium, modelled as a dynamic periodic boundary (Smith et al., 2008).
Furthermore, contact dynamics are also experimentally known to be significant for boundary accumulation behaviours of microswimmers (Kantsler et al., 2013;Bianchi et al., 2017) and to vary extensively with solutes and surfaces (Klein et al., 2003). However, the current study does not consider the detailed surface dynamics in the region very close to the boundary since its scope considers universal hydrodynamic interactions, rather than the boundary behaviours for a specific swimmer, solute, and surface.
Even with a simple short-range repulsion, the details of the repulsive force can alter the swimmer dynamics (Lintuvuori et al., 2016;Ishimoto, 2017), while the contact mechanics reflect the specific biological and physical features of the system under investigation. A further generalisation to be considered in artificial colloidal microswimmers is the impact of sedimentation and gyrotaxis due to swimmer density heterogeneity and density offset from the fluid (Das et al., 2020), as well as the chemical and physical mechanisms that drive the colloidal particle (Uspal et al., 2015). Also, the swirling squirmer characterised by a rotlet-dipole singularity is known to exhibit circular motion near a boundary (Ishimoto and Gaffney, 2013), and the impact of the surface topography on such a swimmer warrants future investigations.
In summary, this investigation has used the nearest neighbour regularised boundary element method (Gallagher and Smith, 2018) to numerically explore the hydrodynamic interactions between a spherical tangential squirmer and a spatially oscillating surface topography with an amplitude that is an order of magnitude less than the squirmer size. In particular, a squirmer was investigated that swam with very simple dynamics close to a flat boundary, relaxing to a stable distance from the wall, and swimming in the direction of its orientation in the horizontal plane. However, even with small amplitude surface topographies, this squirmer's dynamics has depended in a subtle and complex manner on the wavelength of the surface topography. We found that surface topographies could effect limited and local squirmer guidance towards topography troughs, in particular once the squirmer size is of the same order of magnitude as the surface topography wavelength. However, contact dynamics may also be induced at such wavelengths of the surface topography, especially if the initial squirmer orientation to the surface is not along topography crests or troughs. However, surface guided behaviours are robust to other aspects of the surface topography, such as reductions in the amplitude and changes in the shape of the surface topography waves. More generally, the framework enables these predictions to be made forming a basis for in silico experimentation of microorganisms and designing artificial micromachines.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors without undue reservation. The data presented in this article is available from http://dx.doi.org/10.5287/ ora-mzvppgob8.