Hydrodynamics of flagellated microswimmers near free-slip interfaces

TL;DR

This study uses numerical simulations to explore how flagellated microorganisms swim near free-slip interfaces, revealing different behaviors compared to no-slip walls and implications for biofilm formation.

Contribution

It provides new insights into microorganism hydrodynamics near free-slip surfaces, highlighting differences from no-slip boundaries and their effects on swimming trajectories.

Findings

Microorganisms can escape or collide with free-slip surfaces depending on initial conditions.

Swimmers follow a counter-clockwise trajectory near free-slip interfaces.

Hydrodynamics near free-slip surfaces differ significantly from no-slip walls.

Abstract

The hydrodynamics of a flagellated microorganism is investigated when swimming close to a planar free-slip surface by means of numerical solu- tions of the Stokes equations obtained via a Boundary Element Method. Depending on the initial condition, the swimmer can either escape from the free-slip surface or collide with the boundary. Interestingly, the mi- croorganism does not exhibit a stable orbit. Independently of escape or attraction to the interface, close to a free-slip surface, the swimmer fol- lows a counter-clockwise trajectory, in agreement with experimental find- ings, [15]. The hydrodynamics is indeed modified by the free-surface. In fact, when the same swimmer moves close to a no-slip wall, a set of initial conditions exists which result in stable orbits. Moreover when moving close to a free-slip or a no-slip boundary the swimmer assumes a different orientation with respect to its trajectory. Taken together, these results contribute to shed light on the hydrodynamical behaviour of microorgan- isms close to liquid-air interfaces which are relevant for the formation of interfacial biofilms of aerobic bacteria.

Equations

$${\dot{\mathbf{x}}}_J$$

$$\dot{\varphi }$$

$$\dot{\mathbf{q}}$$

$$\mathbf{S}(\mathbf{q})=\left[{\begin{array}[]{rrr}-q_{1}&-q_{2}&-q_{3}\\ q_{0}&-q_{3}&q_{2}\\ q_{3}&q_{0}&-q_{1}\\ -q_{2}&q_{1}&q_{0}\\ \end{array}}\right]\,.$$

$${\mathbf{e}}_1=\left(q_0^2+q_1^2-q_2^2-q_3^2,2(q_1{q}_2+q_0{q}_3),2(q_1{q}_3-q_0{q}_2)\right),$$

$${\mathbf{e}}_2=\left(2(q_1{q}_2-q_0{q}_3),q_0^2-q_1^2+q_2^2-q_3^2,2(q_2{q}_3+q_0{q}_1)\right),$$

$${\mathbf{e}}_3=\left(2(q_1{q}_3+q_0{q}_2),2(q_2{q}_3-q_0{q}_1),q_0^2-q_1^2-q_2^2+q_3^2\right).$$

$$\boldsymbol{\nabla}\cdot\mathbf{u}&=0\,\,,\\ \nabla^{2}\mathbf{u}-\boldsymbol{\nabla}p=0\,\,,$$

$$\int_{H\cup T}{\mathbf{f}dS}&=0\,\,,\\ \int_{H\cup T}{\mathbf{r}\wedge\mathbf{f}dS}=0\,\,.$$

$${\mathbf{e}}_{\mathbf{T}}\cdot {\int }_T\mathbf{r}\wedge \mathbf{f}dS=-{\tau }_M.$$

$$\mathbf{u}(\mathbf{x})&=\mathbf{U}+\mathbf{\Omega}_{H}\wedge\mathbf{r},\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(\mathbf{x}\in H)\,\,,\\ \mathbf{u}(\mathbf{x})=\mathbf{U}+(\mathbf{\Omega}_{H}+\Omega_{T}\mathbf{e}_{T})\wedge\mathbf{r},\,\,\,\,\,(\mathbf{x}\in T)\,\,.$$

$$\mathbf{u}\cdot\mathbf{n}&=0\,\,,\\ (\left(\mathbf{\nabla}\mathbf{u}+(\mathbf{\nabla}\mathbf{u})^{\mathrm{T}}\right)\cdot\mathbf{n})(\mathbf{I}-\mathbf{n}\otimes\mathbf{n})=0\,\,.$$

$$E({\mathbf{x}}_0)u_i({\mathbf{x}}_0)={\int }_{\partial \mathcal{D}}\left[u_j(\mathbf{x})T_{ikj}(\mathbf{x},{\mathbf{x}}_0)n_k(\mathbf{x})-G_{ij}(\mathbf{x},{\mathbf{x}}_0)f_j(\mathbf{x})\right]dS$$

$$\dot{\hat{\Theta }}$$

$$\dot{\hat{h}}$$

$$R=\frac{|\hat{\mathbf{U}}(\hat{\mathbf{\Theta }},\hat{\mathbf{h}})|}{{\hat{\Omega }}_Z-{\hat{\Omega }}_X\tan \hat{\Theta }},$$

$$\mathbf{r}$$

$${\mathbf{r}}^{\prime }$$

$$G_{ij}^{free}=\left(\frac{{\delta }_{ij}}{r}+\frac{r_i{r}_j}{r^3}\right),$$

$$G_{i1}^{FS}=\left(\frac{{\delta }_{i1}}{r}+\frac{r_i{r}_1}{r^3}\right)+\left(\frac{{\delta }_{i1}}{r^{\prime }}+\frac{r_i^{\prime }{r}_1^{\prime }}{r^{\prime 3}}\right),$$

$$G_{i2}^{FS}=\left(\frac{{\delta }_{i2}}{r}+\frac{r_i{r}_2}{r^3}\right)+\left(\frac{{\delta }_{i2}}{r^{\prime }}+\frac{r_i^{\prime }{r}_2^{\prime }}{r^{\prime 3}}\right),$$

$$G_{i3}^{FS}=\left(\frac{{\delta }_{i3}}{r}+\frac{r_i{r}_3}{r^3}\right)-\left(\frac{{\delta }_{i3}}{r^{\prime }}+\frac{r_i^{\prime }{r}_3^{\prime }}{r^{\prime 3}}\right).$$

$$u_i(\mathbf{x})=\frac{1}{8\pi \mu }{G}_{ij}^{FS}(\mathbf{x},\mathbf{X})b_j.$$

$$T_{ik1}^{FS}=-6\frac{r_i{r}_1{r}_k}{r^5}-6\frac{r_i^{\prime }{r}_1^{\prime }{r}_k^{\prime }}{r^{\prime 5}},$$

$$T_{ik2}^{FS}=-6\frac{r_i{r}_2{r}_k}{r^5}-6\frac{r_i^{\prime }{r}_2^{\prime }{r}_k^{\prime }}{r^{\prime }},$$

$$T_{ik3}^{FS}=-6\frac{r_i{r}_3{r}_k}{r^5}+6\frac{r_i^{\prime }{r}_3^{\prime }{r}_k^{\prime }}{r^{\prime }},$$

Full text

1 Introduction

The recent development of micro-fluidic based devices has pushed the scientific community to reconsider classical issues in fluid dynamics such as the hydrodynamical behavior of micro-scale bodies in the creeping regime. For instance, propulsion of microswimmers, see e.g. fauci2006biofluidmechanics ; lauga2009hydrodynamics ; guasto2012fluid ; elgeti2014physics , has been addressed for the production of energy oriented to micro-devices (di2010bacterial, ; sokolov2010swimming, ) and for the self-propulsion of micro-robots (pak2012micropropulsion, ; alouges2013optimally, ). The microswimmer’s behavior in presence of confining surfaces is also crucial for biofilm formation. Most studies (pratt1998genetic, ; rusconi2010laminar, ) focused on the biofilm formation at solid-liquid interfaces. Biofilm formation near air-liquid interfaces is also a problem of significant concern, since the liquid-air biofilm can be advantageous for aerobic microorganisms, providing them access to oxygen (from the air) and nutrients (from the liquid) at the same time (constantin2009bacterial, ; koza2009characterization, ). Several species of flagellated bacteria are known to form biofilms at liquid/air interfaces, e.g. Bacillus subtilis, Bacillus cereus and Pseudomonas fluorescens. For instance, the preparation of the well known Natto from Japanese cooking, consisting of fermented soy beens, involves B. subtilis in significant concentrations (chantawannakul2002characterization, ). In contrast, some strains of B. cereus are known to be harmful to humans and cause foodborne illness (ehling2004bacillus, ). P. fluorescens, instead, is able to contaminate heparinized saline flushes used in cancer therapy (gershman2008multistate, ). Despite its biological implications, the mechanics of flagellated swimmers close to a liquid-air interface is much less explored than the case of a liquid-solid interface.

At micro-scales viscous effects overwhelm inertia leading to the development of apparently counterintuitive swimming strategies as proved by the Scallop theorem (purcell1977life, ; lauga2011life, ). In a nutshell, if the swimmer deforms through a sequence of body configurations which are periodic and time-reversible (reciprocal motion), its average motion is zero. The reciprocal motion can be exploited for locomotion only when non-linear or memory effects are relevant. Examples are the non-Newtonian behavior of the fluid (qiu2014swimming, ; lauga2009life, ; keim2012fluid, ) or motions occurring close to deformable interfaces (trouilloud2008soft, ). Microorganisms developed many different strategies based on nonreciprocal effects to overcome the above restrictions. The Spiroplasma deforms its cytoskeleton by propagating pairs of kinks (trachtenberg2003bacterial, ; yang2009kinematics, ). Other microswimmers exploit their cilia wavy motion (maxey2011biomimetics, ) as done by Paramecium (jana2012paramecium, ).

However, many microorganisms take advantage of single or multiple flagella, such as the Caulobacter crescentus, which has a single (right-handed) helical filament driven by a rotary motor (li2006low, ) and the Escherichia coli that has multiple flagella (berg2004coli, ). The flagellar motor which activates and controls the filament rotation is able to switch between both rotation directions (wang2014switching, ), and, as a first approximation, the torque applied to the filament is constant (berg2004coli, ). In the case of the Escherichia coli the flagella arrange in a bundle, characterized by a left-handed rotation of the motor which confers to bacteria a smooth forward motion (run phase). The flagella is also able to invert the rotation direction in order to let the bacteria at rest and change its direction (tumble phase). Due to their internal structure, the filaments can only assume twelve prescribed shapes (helical polymorphic states) but only one of those is the “normal state”, i.e. it is the most observed one during the run phase (darnton2007force, ; vogel2010force, ). Hence, the flagella can be modeled as a rigid single filament (phan1987boundary, ; shum2010modelling, ) rotating around the bacterial head axis.

The main purpose of this work is an extensive analysis of the hydrodynamical behavior of a flagellated microswimmer close to air-liquid interfaces, with particular attention to free surfaces. For the sake of definiteness, we focus on a simplified geometry modeling the configuration of E. coli whose swimming is probably the most widely studied from both experimental and numerical point of view. Anyhow, in order to explore the effect of different geometrical parameters, several modified configurations are also discussed.

Flagellated swimmers close to a free-surface have been already addressed by simplified models, see e.g.  crowdy2011two ; di2011swimming ; lopez2014dynamics where two-dimensionality, resistive force theory, and multipole expansions techniques were exploited, respectively. To the best of our knowledge, the present work addresses the first fully three dimensional numerical simulation of a swimmer in presence of a liquid-air interface. Given the linearity of the Stokes equations which are the appropriate model for creeping flows, the numerical approach exploits the Boundary Element Method (BEM). The BEM can easily handle the complex geometry of the microswimmer and account for different boundary conditions. Moreover, a flat and infinitely extended free surface and/or solid wall (when needed for comparisons) can be easily included by considering the appropriate Green’s function, thus avoiding any undesired effect of the numerical truncation of the domain.

There is a wide body of literature available dealing with the motion of microorganisms in free space or close to solid boundaries. Since the first studies by taylor1951analysis , the attention was principally focused on microorganisms whose flagellum was modeled as an infinite cylindrical filament in an unbounded fluid endowed with small amplitude wavy motion. In fact, it was initially believed that the flagella were only moved by wave propagation. However berg1973bacteria showed that bacteria could also rotate their flagella in a corkscrew-like motion, moving the flagellar bundle as a single filament. The flagellum has a large aspect ratio, with length exceeding thickness by even more than two orders of magnitude (lauga2009hydrodynamics, ). This particular geometry pushed the adoption of a Slender Body Theory (SBT) that has been exploited in several studies to evaluate the translational velocity and/or the torque applied by the swimmer (hancock1953self, ). Successively higdon1979hydrodynamic ; higdon1979hydrodynamics transformed the Stokes equations into a system of singular integral equations accounting for the swimmer translational and angular velocities. He added to the SBT a variable strength of the singularity along the flagellum centerline, thus modelling different centerline geometries. In this case both the planar sinusoidal motion and the rotation about the body axis were amenable to modelling. By an appropriate system of images the SBT could also account for a spherical cell body and the presence of a wall. However too many restrictions still confined the SBT application to extremely simplified configurations.

In recent years the increase of computational resources, led to a more extensive use of the BEM which overcomes many drawbacks of the SBT in dealing with microswimmers both in free space and confined conditions. phan1987boundary used the BEM to study the motion of a microorganism in free space. Successively, ramia1993role applied the BEM to the interaction between the swimmer and a solid wall, showing that, when swimming close to a solid wall, the swimmer exhibited a circular motion. The BEM was also used to study the interactions between two neighboring flagellated microswimmers highlighting the possibility of a coordination between their flagellar motion in order to maximize their velocities (ramia1993role, ).

More recently, lauga2006swimming experimentally investigated the motion of an E. coli near a solid wall and found, as predicted by ramia1993role , a circular clockwise motion. In the paper the authors also provided a simple theoretical model which was able to explain their experimental observations. Interestingly, the same model suggested that the bacteria should reverse its rotation when swimming in proximity of a free surface. The same authors also showed that the swimmer is attracted by a solid wall. This behavior was deeply investigated a few years later by giacche2010hydrodynamic in a numerical work which highlighted how the bacteria could move at a stable distance from a wall. At the same time shum2010modelling investigated the motion of a microswimmer close to solid walls by considering many geometrical configurations and relating the shape of the body to the propulsion efficiency and to the possibility of achieving different motions, i.e. to be attracted by the wall, to escape from the wall or to reach stable circular orbits at a given wall-normal distance. The same authors (shum2012effects, ) also focused on the flexible hooks linking each flagellum to the cell. They studied the modifications in the microswimmer trajectories when changing the hook rigidities. They found that, within an intermediate range of rigidities, the swimmer behavior doesn’t change too much with respect to the simpler model of rigid hook. The same work highlighted how, for particular values of relative hook stiffness, there is a transient phase of periodic motion with constant average distance from the wall, leading to boundary accumulation. The tendency of swimming microorganisms to accumulate near solid walls (li2009accumulation, ), with particular emphasis on collisions with the surface and rotational Brownian motion, has also been investigated and linked to the swimming speed and the cell size (li2011accumulation, ) of the microswimmer. Recent works extended the investigation about the motion of microswimmers to more complex surfaces: a clean fluid-fluid interface, a slipping rigid wall, and a fluid interface covered by surfactants (lopez2014dynamics, ). The authors used an asymptotic, far field approximation to represent the actual swimmer, retaining information about velocities and rotations. The case of two fluid interfaces, in the limit of vanishing viscosity of one of them, allowed to describe a free surface. In such conditions, the swimmer exhibited counter-clockwise motion. This confirmed the results anticipated by lauga2006swimming . Also di2011swimming supported by means of experimental observations the theoretical prediction made through a simplified model based on the method of images and the resistive force theory. Even if neglecting all the dynamics of the swimmer outside the interface plane, this work provided a simple explanation for a counter-clockwise motion over a perfectly slipping surface, showing a good agreement between the results of the simplified model and the experimental observation. Briefly,  lauga2006swimming explain that, when swimming above a no slip surface, a positive rotation rate of the swimmer head around its longitudinal axis produces a lateral force which is opposite to the one induced by the negative tail rotation rate. As a consequence, a net torque normal to the wall is exerted on the body such that the swimmer follows a clockwise (CW) trajectory. Otherwise, when swimming close to a liquid-air interface, di2011swimming devise a simplified model based on resistive force theory endowed with suitable symmetries to satisfy the free slip condition at the interface. In this model, the swimmer moves under the effect of the velocity generated by its mirror image below the interface. The counter-rotating image head produces on the swimmer head a lateral velocity which is opposite to the force that is exerted in the case of a solid wall. Such relative velocity gives rise to a corresponding viscous force in the same direction. Since the same reasoning applies to the counter rotating tail, the overall torque on the microswimmer is also opposite to the one experienced on a no-slip wall, hence a net CCW motion is produced. In literature a clockwise motion has been observed (lemelle2010counterclockwise, ) also in presence of a free surface. Based on experimental observations, morse2013molecular attributed such a behavior to the molecular adsorption (due to the presence of biological material in the growth medium) altering the rheological properties of the air/water interface, thus determining the swimming pattern of nearby cells.

In principle, the mechanical causes that can affect the swimming direction are: $i)$ the rotation direction of the flagellum bundle; $ii)$ the effective boundary condition at the planar surface. This study focuses on the point $ii)$ through the numerical simulation of the motion of a E. coli-like microswimmer close to free-slip and no-slip surfaces, assuming a standard left-handed arrangement for the flagellum bundle. The aim is investigating the behavior of the swimmer by addressing in full detail the complete three dimensional nature of the hydrodynamical interaction between the swimmer and the surface. It is worth noting that free-slip and no-slip are the limits of the more general Navier boundary condition that connects the velocity at the liquid boundary with the tangential stress (bazant2008tensorial, ). For liquid water moving on a solid surface, the actual slip is negligible at micro-scale also for hydrophobic coatings (chinappi2010intrinsic, ; sega2013regularization, ) and only the presence of vapor bubbles trapped in the surface asperities (superhydrophobic surfaces) leads to significant slippage (gentili2014pressure, ; bolognesi2013novel, ), potentially altering the motion of particles close to the surface  (pimponi2014mobility, ; nizkaya2014flows, ). Hence, the no-slip boundary is, for the present purposes, a reliable model of a rigid wall. On the other hand, the proper boundary conditions at the liquid-air interface are impermeability and continuity of the tangential stress components, that reduces to free-slip (zero tangential stress) since air density is order of magnitude smaller than the liquid one.

The paper is organized as follows: the geometrical model of the swimmer and the BEM is addressed in § 2. § 3 reports the salient results concerning the motion of the swimmer in presence of a free surface. Finally, we discuss the major findings and point out the main conclusions (§ 4), giving a perspective view of future work.

2 Boundary integral formulation for the microswimmer

2.1 Swimmer model

This section concerns the modelling of a microswimmer inspired to Escherichia coli. This bacteria has been deeply investigated in the literature and a lot of information about its geometry and propulsion mechanism is available. E. coli has a cell length which varies between $1.6$ and $3.9\,\mu m$, the cell width ranges between $0.9$ and $1.7\,\mu m$, with a resulting cell volume from $1.5$ to $4.4\,\mu m^{3}$ (volkmer2011condition, ) and an average length of flagella of about $7\mu m$ (lauga2006swimming, ).

Figure 1 sketches the simplified geometry comprising the ellipsoidal axisymmetric cell and the corkscrew tail. Hereafter the equivalent radius of the ellipsoidal cell $\bar{a}$ (the radius of the sphere having the same volume) will be used as reference length-scale, i.e. the dimensionless cell volume is $V=3V^{\prime}/{4\pi\bar{a}}^{3}=1$, where the prime identifies dimensional quantities. The cell aspect ratio is $AR=a_{1}/a_{2}=2$, being $2a_{1}$ and $2a_{2}$ the longitudinal and the transversal (dimensionless) axis. Following shum2010modelling , the tail bundle is modelled as a single helix with radius $a_{T}=0.05$. The tail rotates around its axis $\mathbf{e_{T}}$ (see figure 1). The dimensionless axial length of the helix is $L$, $A$ denotes the helix amplitude and $\lambda$ its wavelength. In the present work we selected typical values for the tail length and cell axes, namely $L=7$, $a_{1}=1.6$ and $a_{2}=0.8$, while the amplitude $A$ and the wavelength $\lambda$ of the tail are systematically varied. With the above choices, the swimmer is rescaled into an actual E. coli by assuming $\bar{a}=1\,\,\mu m$.

It is useful to introduce a body reference frame with orthonormal base vectors $\{\mathbf{e_{1}},\mathbf{e_{2}},\mathbf{e_{3}}\}$ where $\mathbf{e_{1}}=-\mathbf{e}_{T}$ is longitudinal (see figure 1). The subscript $H$ identifies the head and the subscript $T$ refers to the tail, which can rotate with respect to the head about its axis $\mathbf{e_{T}}$. During the tail rotation around $\mathbf{e}_{T}$, each point of the rigid tail describes a circle in the plane normal to $\mathbf{e}_{T}$, being $\phi$ the corresponding rotation angle (or flagellum phase). The time derivative $\dot{\phi}(t)$ is the tail rotational velocity $\Omega_{T}$. The swimmer position is identified in a fixed reference frame with base $\{\mathbf{X},\mathbf{Y},\mathbf{Z}\}$ by the three coordinates of the cell-to-tail junction point $\mathbf{x}_{J}$, see figure 2. The swimmer translates with dimensionless velocity $\mathbf{U}=\mathbf{U^{\prime}}/v$ and has angular velocity $\mathbf{\Omega}_{H}=\mathbf{\Omega^{\prime}}_{H}\bar{a}/v$, where $v\simeq 20\mu m/s$ is a typical swimming velocity taken as reference quantity. The $\{\mathbf{X},\mathbf{Y}\}$ coordinate plane is taken to coincide with the planar boundary, either the solid wall or the free surface, see figure 2.

The kinematics of the swimmer is described by seven degrees of freedom, the junction position $\mathbf{x}_{J}(t)$, the flagellum phase $\phi(t)$ and the three parameters defining the rotation of the head. We stress that the discretely evolved rotation matrix should belong to the proper matrix subspace, namely the $\rm SO(3)$ subgroup. This is enforced through a description in terms of quaternions with unit norm, see e.g. diebel2006representing . It follows,

[TABLE]

where $\mathbf{q}=(q_{0},q_{1},q_{2},q_{3})$ is the quaternion vector with $|\mathbf{q}|=1$. When $\mathbf{\Omega}_{H}(t)$ is expressed in the $\{\mathbf{X},\mathbf{Y},\mathbf{Z}\}$ base, $\mathbf{S}$ reads

[TABLE]

Given the quaternion $\mathbf{q}(t)$, the components of the body frame unit vectors can be reconstructed as

[TABLE]

Time integration of equations (1),(2), (3), allows to track the swimmer trajectory once $\mathbf{U}$, $\mathbf{\Omega}_{H}$ and $\Omega_{T}$ are determined by the solution of the hydrodynamical interactions between the swimmer and the fluid, as discussed in the next section. Since the numerical integration error in the equations for the quaternion may affect its norm, the quaternion is normalized to unit length at each time step. It easily shown that the accuracy of the integration scheme is exactly preserved by this procedure.

2.2 Hydrodynamic model

Given the characteristic length scale $\bar{a}$ and swimming velocity $v$ of the microorganism, the typical Reynolds number in water is $Re=\rho\bar{a}v/\mu\simeq 2*10^{-5}$ where $\mu$ and $\rho$ are the water viscosity and density, respectively. Hence the inertial terms are negligible and, for Newtonian fluids, the flow in the domain ${\cal D}$ is described by the Stokes equations,

[TABLE]

where $\mathbf{u}=\mathbf{u^{\prime}}/v$ and $p=p^{\prime}\bar{a}/(\mu v)$.

The flow velocity $\mathbf{u}$ is forced by the tension $\mathbf{f}$ at the swimmer surface acting on the fluid. The swimmer propulsion is due to the internal (constant) torque $\tau_{M}=\tau^{\prime}_{M}/(\mu v\bar{a}^{2})$ exchanged between head and tail. Being a free body, the total force and torque on the swimmer (head and tail) is zero,

[TABLE]

where $\mathbf{r}=\mathbf{x}-\mathbf{x}_{J}$. The torque $\tau_{M}$ exerted on the tail is balanced by the torque produced by the fluid stresses $\mathbf{f}$ on the tail boundary, namely

[TABLE]

The system  (6),  (7) and  (8) needs proper boundary conditions at the swimmer surface and external boundaries, i.e. wall or free surface. On the microswimmer surface the no-slip condition yields

[TABLE]

No-slip boundary condition is also used to model the solid wall. For the free surface vanishing normal fluid velocity (impermeability) and zero tangential stresses (free slip condition) are the appropriate prescriptions,

[TABLE]

2.3 Boundary Element Method

The solution of the Stokes equations (6) can be recast in integral form

[TABLE]

where $E(\mathbf{x_{0}})=1$ for points belonging to the interior of the domain. $G_{ij}(\mathbf{x},\mathbf{x}_{0})$ is the free-space Green function, i.e. the $i-$th component of the velocity at $\mathbf{x}$ induced by a Dirac $\delta$ like force at $\mathbf{x}_{0}$ acting in direction $j$, and $T_{ikj}(\mathbf{x},\mathbf{x_{0}})$ the associated stress tensor (happel1983low, ). Equation (11) expresses the $i$-th velocity component at the collocation point $\mathbf{x}_{0}$ in the fluid domain ${\cal D}$. It requires the knowledge of the velocity $u_{j}(\mathbf{x})$ and the stresses $f_{j}(\mathbf{x})$ at the boundary $\partial{\cal D}$. The integral representation (11) can be turned into a boundary integral equation in the limit $\mathbf{x}_{0}\to\partial{\cal D}$. The resulting expression can be written in the same form (11) now with $E(\mathbf{x}_{0})=1/2$ and the integral understood in the Cauchy principal value sense. We emphasize that the boundary $\partial{\cal D}$ consists of two disjointed parts: the swimmer surface and the planar interface. The integral can be restricted to the swimmer’s surface by exploiting the symmetries associated with the boundary condition at the planar interface. This is tantamount to the use of the appropriate Green’s function, see appendix A, for the free-slip case and blake1971note for the no-slip case.

Due to the boundary condition (9) the fluid velocity at the swimmer surface is expressed in term of seven unknowns, namely, the velocity of the junction point $\mathbf{U}$, the angular velocity $\mathbf{\Omega}_{H}$, and the tail rotation velocity $\Omega_{T}$, which together with the tension $\mathbf{f}$ at the swimmer surface complete the set of unknowns. The system of equations consists of the vector boundary integral equation (11), the global force and torque balance (7), and the torque balance for the swimmer tail (8). It can be shown that the solution exists and is unique (ladyzhenskaya1969mathematical, ).

The above system is discretized by means of $N$ curved 6-point elements (typically $N=3518$ total elements, with $N_{H}=512$ panels on the swimmer head), as sketched in figure 2, with piecewise constant shape functions (pozrikidis1992boundary, ). By selecting the center of each element as a collocation point $\mathbf{x_{0}}$, the boundary integral equation (11) with $E(\mathbf{x_{0}})=1/2$ is recast in a system of $3N$ scalar equations for the stresses. As will be shown in § 3 the swimmer may approach the interface. In this case some of the panels belonging to the swimmer and to its image may get very close. This potentially spoils the accuracy of the integrals providing the influence coefficients appearing in the discrete equations. Care is taken in the simulations discussed below to have the diameter of the typical panel smaller than the distance between the body and its image thereby preventing the undesired loss of accuracy. Concerning the solution of the algebraic system, inverting a matrix coming from the discretization of boundary integral equations may require some care, since the algebraic system may be prone to severe ill-conditioning. An example is provided by the piecewise constant approximation of a first kind Fredholm operator (hsiao1973solution, ). In the context of the Stokes equations, such operator arises when solving for the stresses given the velocity. However, the problem of the swimmer addressed in the present paper does not involve a pure first kind Fredholm operator, since the equations are augmented by the free-body constrains (zero net forces and torques) and by the equation enforcing the internal torque, with unknowns the stresses on the body plus the rigid body velocities and the relative tail-body rotation rate. The spectral properties of the ensuing matrix differ substantially from those of a pure first kind Fredholm equation. In fact, the analysis based on the singular value decomposition, (golub2012matrix, ), excludes matrix ill-conditioning. As a consequence, the complete system of $3N+7$ algebraic equations for the $3N+7$ unknowns collectively denoted $\boldsymbol{\chi}$, $\mathbf{A}\cdot\boldsymbol{\chi}=\mathbf{b}$, can be solved by standard techniques, like the Gauss-Jacobi method here implemented in an in-house MPI parallel code. Once the solution $\boldsymbol{\chi}$ is obtained, the integration of the kinematic equations (1), (2),  (3), performed via a third order low-storage Runge-Kutta method, allows to track the swimmer trajectory.

3 Swimming close to a free surface: results and discussion

3.1 Phase plane analysis and trajectories for a no-slip wall

The interaction between the swimmer and a homogeneous planar surface is completely determined by three parameters: the distance $h$ between the junction and the surface, the flagellum rotation phase angle $\phi$ and the pitch angle $\Theta$ between the swimmer longitudinal axis $\mathbf{e_{1}}$ and the surface plane, see figures 1 and 2. $\phi$ is a fast variable, thus, following shum2010modelling , the kinematics of the system can be described in term of $\phi-$averaged quantities, in the following denoted by a circumflex. The kinematics in the reduced $(\hat{\Theta},\hat{h})$ space is ruled by

[TABLE]

where $\hat{\Omega}_{2}(\hat{\Theta},\hat{h})=\hat{\mathbf{\Omega}}_{H}\cdot{\mathbf{e}_{2}}$ and $\hat{U}_{Z}(\hat{\Theta},\hat{h})=\hat{\mathbf{U}}\cdot{\mathbf{Z}}$. Figure 3(a) shows the phase plane for a microswimmer close to a solid no-slip surface. Three reduced trajectories, $I$, $II$ and $III$ are traced. The trajectory $III$ hits the wall while trajectories $I$ and $II$, although starting from different initial configurations, both converge to the same attractor. The attractor is an asymptotically stable equilibrium point (cencini2009chaos, ), ($\hat{\Omega}_{2}$ and $\hat{U}_{Z}$ are zero), i.e. nearby trajectories converge to the equilibrium point.

The $\phi-$averaged trajectory in 3D space corresponding to $I$ is obtained by integrating the $\phi-$averaged version of the kinematic equations (1), (3) and is reported in figure 3(b). Apparently, after a transient, the swimmer stabilizes on a circular clockwise (CW) orbit which corresponds to the stable point in the reduced space $(\hat{\Theta},\hat{h})$. Following shum2010modelling , the curvature radius of the stable orbit is given by

[TABLE]

which yields $R=-22.1$. As expected this result matches within numerical accuracy its direct measure $R=-21.9$ obtained from the 3D trajectory, see figure 3(b). Note that equation (14) differs in sign from the original equation given in shum2010modelling due to a different choice of the body reference frame. It is also worth noting that equation (14) provides the exact curvature radius only for the stable orbit. Nevertheless, it also gives a reasonable approximation in the other conditions here explored, as discussed below. In figure 3(a) the color map refers to $R$ as estimated from eq. (14). The color code in figure 3(b) corresponds to the local curvature radius of the trajectory projection in the $xy$-plane. This information is also reported along the curve $I$ in the phase plane, panel 3(a), through the color of the open circles superimposed to the trajectory $I$. From the data it is apparent that the color differences between the circles and the color map in the background can hardly be appreciated.

3.2 Swimming in presence of a free-slip interface

The phase plane analysis has been performed for a microswimmer moving close to a free-slip interface, as illustrated in figures 3(c) and 3(d). The first significant result is the sign of the curvature radius, now positive, i.e. the swimmer exhibits a counter-clockwise (CCW) motion in contrast to what observed for the no-slip wall, compare the panels (b) and (d) in figure 3. The curvature radius is smaller with respect to the no-slip case at corresponding $(\hat{\Theta},\hat{h})$. The curvature radius measured in our simulation favorably compares with experimental observation of di2011swimming where a CCW motion with a curvature radius of $\simeq 10\mu m$ is described. Indeed, assuming $\bar{a}=1\mu m$, appropriate for an E.coli, the typical radius of curvature found from the present numerics, $R\simeq 7$, gives $R^{\prime}=R\bar{a}\simeq 7\mu m$.

The second significant result is that no stable trajectory exists when the swimmer moves close to a free surface. In fact, a swimmer similar to an E.coli moving almost parallel to the free-surface (low $\hat{\Theta}$) is always attracted to it. This result confirms the indication of the 2D model proposed by crowdy2011two for undeformable interfaces. This conclusion also partially agrees with the results illustrated by lopez2014dynamics , where the authors, basing on a multipole expansion technique, highlighted that the average surface-normal velocity $\hat{U}_{z}$ of the swimmer is always negative for $\hat{\Theta}=0$. Our work adds that the surface is no more able to attract the swimmer if it swims above a certain height $h$, even with zero tilt angle.

In figure  4 we extend our analysis to a wider range of initial conditions, for both no-slip and free-slip surfaces. In both cases, at high positive $\hat{\Theta}$ the microswimmer escapes from the surface while, at high negative $\hat{\Theta}$, it eventually hits the boundary. On the no-slip wall three cases are possible depending on the initial conditions, namely the swimmer achieves a stable orbit, it escapes from the wall or it collides with the wall. In contrast, for a free-slip surface only two possibilities exist, namely, escaping from the surface or being attracted towards the interface. The maps also highlight the presence of a further equilibrium point $(\hat{\Theta}_{u},\hat{h}_{u})$ where the time derivatives of both ${\hat{h}}$ and $\hat{\Theta}$ vanish. However, as shown by the behavior of the streamlines in the neighbor of $(\hat{\Theta}_{u},\hat{h}_{u})$, this point is unstable for both no-slip and free-slip cases, i.e. nearby trajectories escape from the equilibrium point.

In order to extend the analysis, the dynamics of swimmers with different geometrical characteristics was also examined, focusing in particular on the shape of the head and the tail (axial) length. Figure 5 shows the phase plane $\hat{\Theta}-\hat{h}$ when the tail length $L\in(3,5,10,15)$ and the head aspect ratio $AR=a_{1}/a_{2}\in(1,3,4,5)$. When changing the aspect ratio, the swimmer does not substantially modify its behavior. In particular, the orientation as described by the pitch angle $\hat{\Theta}$ does not qualitatively change and the regions of phase space where the swimmer is attracted to the free-surface or where it escapes away from it remain very similar. The general trend is an increase of the curvature of the trajectory for given $\left({\hat{\Theta}},\,{\hat{h}}\right)$ as $AR$ increases. On the other hand, when changing the relative tail-to-head length, the swimmer hits the wall with different pitch angles $\hat{\Theta}$. Decreasing the tail length, the region of phase plane where the swimmer escapes becomes larger while the region where it is attracted to the free-surface shrinks. The curvature of the trajectory, for given $\left({\hat{\Theta}},\,{\hat{h}}\right)$, is found to increase with the tail length.

3.3 Drift along the trajectory

Typically the longitudinal axis of the microswimmer is misaligned with the local velocity evaluated at the junction point. This behavior was experimentally observed for E.Coli moving close to a liquid-air interface by di2011swimming . In order to compare with experiments we define the drift angle $\alpha$ as the angle between the ${X,Y}$ projections of the body frame unit vector $\mathbf{e}_{1}$ and the swimmer velocity $\mathbf{U}$. The sign of $\alpha$ is taken to be positive when the swimmer points outward with respect to the trajectory, see figure 6a.

In figure 6b and c, the $\phi-$averaged drift angle is reported for both no-slip and free-slip cases. In the free-slip case the head points outside the trajectory and $\hat{\alpha}$ increases as the swimmer approaches the wall. The spanned range of values, $\hat{\alpha}\in(10^{\circ},30^{\circ})$, is in good agreement with the experimental observation (di2011swimming, ). In contrast, in the no-slip case, $\hat{\alpha}$ is slightly negative, $\hat{\alpha}\simeq-2^{\circ}$, meaning that the swimmer is almost aligned with the trajectory of the reference point $\mathbf{x}_{J}$.

Here, to complete the discussion, the effect of modifying the tail geometry is briefly addressed. Indeed, among the large number of different parameters defining the flagellated geometry, probably the most uncertain ones concern the tail, specifically amplitude of the helix, $A$, and number of turns, $N_{\lambda}$. Figure 7 reports the radius of curvature $R$ and the drift angle $\hat{\alpha}$ for three different tails. The comparison shows that the results discussed so far are generic.

4 Conclusions

The motion of a flagellated microswimmer close to a boundary, either a solid wall or a free surface, is relevant to several applications spanning from micro-robots to biology and medicine, as concerning in particular biofilm formation. Its dynamics can in principle be affected by several physical phenomena occurring at the microscale and, for biological applications, by the behavior of the microorganism. The chemical and physical nature of the interface may play a role in the interaction between the swimmer and the surface, e.g. surface charges and chemicals adsorbed at the interface may have a significant influence. Concerning the specific case of E. coli, taken as representative of most flagellated, recent experimental studies reported a characteristic motion of the swimmer in the two extreme cases of a solid wall and a free surface. In the first case, the experimental observation consistently show that the microswimmer typically moves in circulatory orbits oriented in clockwise direction (CW) lauga2006swimming . In contrast, there is evidence that the orientation of the trajectory is reversed (counterclockwise, CCW) when swimming occurs near a free surface di2011swimming . However, the effect of the free surface is less neatly defined and, sometimes, CW motion is reported, probably due to the presence of contaminants adsorbed at the interface (lemelle2010counterclockwise, ). Theoretical models proposed to explain this behavior typically exploit several form of approximation, e.g. multipole expansion, that may become partially inaccurate when the microswimmer gets very close to the interface.

The present paper provides a complete description of the motion, considering a reasonably realistic geometry of the flagellated in presence of a free-surface modeled as a rigid, free-slip plane. The resulting model removes any concurrent effect, retaining a full hydrodynamics description. The results for the free-slip boundary were compared with those already analyzed in shum2010modelling for the no-slip case. The data clearly indicate that the motion close to a liquid-air interface is CCW, in agreement with the experimental data (di2011swimming, ) and with the theoretical results obtained by using multipole expansions (lopez2014dynamics, ) and resistive force theory (di2011swimming, ).

Other available experimental information, namely the orientation of the bacteria with respect to its trajectory (di2011swimming, ), is satisfactorily reproduced by the present simulations, confirming that the head of the swimmer points outward. In contrast, the bacteria is roughly aligned with its trajectory close to a no-slip wall. In principle, this observable can be used together with the rotation direction to interpret experimental results on the interaction of a microswimmer with an interface. A characteristic aspect of the motion near a solid surface is the occurrence of stable orbits (giacche2010hydrodynamic, ; shum2010modelling, ). To the contrary, no stable orbit has been presently found on a free-slip interface.

In conclusion, the boundary conditions are confirmed to deeply influence the hydrodynamical behavior of the swimmer. This consideration paves the way to suitably textured surfaces which, properly engineered to stably support a super-hydrophobic state, can be exploited to passively control the microswimmer motion. Indeed the ability of superhydrophobic surfaces to alter the flow has been recently used for passive particle separation (asmolov2015principles, ). Under this respect fully resolved hydrodynamic simulations able to model complex physical surfaces, like the BEM adopted here, could provide the required fundamental knowledge to extend passive control strategies to active suspensions.

The authors acknowledge the CINECA Iscra C Award (IscrC-BSM-LAI) for the availability of HPC resources and the ERC Grant No. [339446] BIC: Bubbles from Inception to Collapse.

Appendix A Green functions for a flat free-slip interface

The Green’s function for a solid adherent wall is known since the work of blake1971note . The corresponding fundamental solution for a free-slip interface does not seem to be mentioned in the current literature. As shown here, it can be obtained by means of the method of images where the generic localized force is reflected with respect to the plane in such a way to satisfy the impermeability and the zero tangential stress condition. For convenience the index notation is adopted here with $\mathbf{x}=(x_{1},x_{2},x_{3})$ the coordinates of a point in space. The image system is reported in figure 8 (left panel) where the point force $\mathbf{b}=(b_{1},b_{2},b_{3})$ placed at $\mathbf{X}=(X_{1},X_{2},X_{3})$ is reflected in $\mathbf{X^{\prime}}=(X_{1},X_{2},-X_{3})$ as $\mathbf{b^{\prime}}=(b_{1},b_{2},-b_{3})$. The relative positions

[TABLE]

connect the generic point $\mathbf{x}$ in the fluid domain with the two singularities in $\mathbf{X}$ and $\mathbf{X}^{\prime}$, respectively. Taking into account the well known expression for the free space Green function (happel1983low, )

[TABLE]

the superposition of the two singularities yields

[TABLE]

The i-th velocity component at $\mathbf{x}$ is expressed in terms of the force intensity $\mathbf{b}$ at $\mathbf{X}$ as

[TABLE]

It can be immediately checked that the ensuing velocity field satisfies impermeability and free-slip conditions at $x_{3}=0$. Indeed, as concerning the wall normal velocity, impermeability straighforwardly follows from (21). Let’s consider now the zero tangential stress condition. From the stress tensor associated with (21)

[TABLE]

one immediately realizes that the tangent component of the tension at $x_{3}=0$ vanishes identically, i.e. $T_{i3j}^{FS}b_{j}=0$ for $i=1,2$. An illustration is provided in the two rightmost panels of figure 8 where the streamlines associated with free-slip Green function due to a point force respectively parallel $\mathbf{b}=(1,0,0)$ and normal $\mathbf{b}=(0,0,1)$ to the wall placed at $X=(0,0,1)$ are reported.