Model microswimmers in channels with varying cross section

TL;DR

This paper develops an analytical framework to understand how microswimmers behave in channels with varying cross sections, revealing how geometry and swimming mechanism influence accumulation and rectification.

Contribution

It extends the Fick--Jacobs approach to microswimmers in corrugated channels, providing new insights into their distribution and rectification based on swimmer type and channel geometry.

Findings

Microswimmer accumulation at walls depends on swimming mechanism and channel shape.

Effective equations describe swimmer distribution along channels with weak corrugation.

Asymmetric channels induce rectification of microswimmers, influenced by swimmer type.

Abstract

We study different types of microswimmers moving in channels with varying cross section and thereby interacting hydrodynamically with the channel walls. Starting from the Smoluchowski equation for a dilute suspension, for which interactions among swimmers can be neglected, we derive analytic expressions for the lateral probability distribution between plane channel walls. For weakly corrugated channels we extend the Fick--Jacobs approach to microswimmers and thereby derive an effective equation for the probability distribution along the channel axis. Two regimes arise dominated either by entropic forces due to the geometrical confinement or by the active motion. In particular, our results show that the accumulation of microswimmers at channel walls is sensitive to both, the underlying swimming mechanism and the geometry of the channels. Finally, for asymmetric channel corrugation our model predicts a rectification of microswimmers along the channel, the strength and direction of which strongly depends on the swimmer type.

Equations

$$\frac{\partial }{\partial t}\rho (\mathbf{r},\theta )=-\vec{\nabla }\cdot \mathbf{J}(\mathbf{r},\theta )$$

$$\mathbf{J}=\left(\begin{array}[]{c}\mathbf{J}_{\mathbf{r}}\\ J_{\theta}\end{array}\right)\qquad\mathrm{and}\qquad\vec{\nabla}=\left(\begin{array}[]{c}\vec{\nabla}_{\mathbf{r}}\\ \frac{\partial}{\partial\theta}\end{array}\right)$$

$${\mathbf{J}}_{\mathbf{r}}(\mathbf{r},\theta )$$

$$J_{\theta }(\mathbf{r},\theta )$$

$$C(\mathbf{r})=C(x,z)=\{\begin{array}{ll}0& |z|\le h(x)-R+\mathcal{O}(|{\partial }_{x}h(x){|}^2)\\ \infty & |z|>h(x)-R+\mathcal{O}(|{\partial }_{x}h(x){|}^2),\end{array}$$

$$\vec{\nabla }\cdot \mathbf{J}={\vec{\nabla }}_{\mathbf{r}}\cdot {\mathbf{J}}_{\mathbf{r}}+\frac{\partial }{\partial \theta }{J}_{\theta }=0,$$

$${\nabla }_{\mathbf{r}}^2\rho +{\vec{\nabla }}_{\mathbf{r}}\cdot [\rho {\vec{\nabla }}_{\mathbf{r}}\beta C]-{\mathrm{Pe}}_t{\vec{\nabla }}_{\mathbf{r}}\cdot \left[\left(\mathbf{n}+\mathbf{v}\right)\rho \right]+\frac{{\mathrm{Pe}}_t}{{\mathrm{Pe}}_r}\frac{\partial }{\partial \theta }\left[\left(\frac{\partial }{\partial \theta }-\frac{{\omega }_0}{D_r}\omega \right)\rho \right]=0,$$

$${\mathbf{n}}_h(x)\cdot {{\mathbf{J}}_{\mathbf{r}}|}_{\pm \frac{h(x)-R}{d}}=0$$

$$\rho (\mathbf{r},\theta )=\rho (\mathbf{r}+L{\mathbf{e}}_x,\theta )$$

$$\int d\mathbf{r}d\theta \rho (\mathbf{r},\theta )=1$$

$$J_{\theta }{|}_{\theta =0}=J_{\theta }{|}_{\theta =2\pi }$$

$${\int }_0^{2\pi }d\theta {\int }_{\frac{-h(x)+R}{d}}^{\frac{h(x)-R}{d}}dz{J}_{\theta }(\mathbf{r},\theta )=0$$

$$J_z(z=\pm 1,\theta )$$

$${\nabla }_{\mathbf{r}}^2\rho +{\nabla }_{\mathbf{r}}\cdot \left[\rho {\nabla }_{\mathbf{r}}\beta C\right]-{\mathrm{Pe}}_t{\nabla }_{\mathbf{r}}\cdot \left[\left(\mathbf{n}+\mathbf{v}\right)\rho \right]=0$$

$$\frac{\partial }{\partial z}\rho (z,\theta )-\left({\mathrm{Pe}}_t[\cos \theta +v_z(z,\theta )]+\frac{\partial }{\partial z}\beta C(z)\right)\rho (z,\theta )=-J_{z,0}(\theta )=0.$$

$${\rho }_0(z,\theta )={\lambda }_0(\theta )e^{{\mathrm{Pe}}_t\left[z\cos \theta +\int v_z(z,\theta )dz\right]+\beta C(z)}={\lambda }_0(\theta )e^{-{\mathrm{Pe}}_{t}W(z,\theta )}$$

$$W(z,\theta )=\{\begin{array}{ll}-z\cos \theta -\int v_z(z,\theta )dz& |z|\le 1\\ \infty & |z|>1\end{array}$$

$$\int dz{J}_{\theta }(z,\theta )=0$$

$$\frac{1}{D_r}{\int }_{-\infty }^{\infty }{J}_{\theta }(z,\theta )dz=-{\lambda }_0\frac{\partial }{\partial \theta }{\int }_{-\infty }^{\infty }{\hat{\rho }}_{0}dz-{\int }_{-\infty }^{\infty }{\hat{\rho }}_{0}dz\frac{\partial }{\partial \theta }{\lambda }_0+\frac{{\omega }_0}{D_r}{\lambda }_0{\int }_{-\infty }^{\infty }\omega {\hat{\rho }}_{0}dz=0$$

$${\lambda }_0(\theta )={\bar{\rho }}_0\frac{1}{{\int }_{-\infty }^{\infty }{\hat{\rho }}_{0}dz}\exp \left[\frac{{\omega }_0}{D_r}\int d\theta \frac{{\int }_{-\infty }^{\infty }\omega (z,\theta ){\hat{\rho }}_0(z,\theta )dz}{{\int }_{-\infty }^{\infty }{\hat{\rho }}_0(z,\theta )dz}\right]$$

$${\rho }_0(z,\theta )={\bar{\rho }}_0\frac{e^{-{\mathrm{Pe}}_{t}W(z,\theta )}}{{\int }_{-\infty }^{\infty }{e}^{-{\mathrm{Pe}}_{t}W(z,\theta )}dz}\exp \left[\frac{4}{3}\frac{R}{d}{\mathrm{Pe}}_t\int d\theta \frac{{\int }_{-\infty }^{\infty }\omega (z,\theta )e^{-{\mathrm{Pe}}_{t}W(z,\theta )}dz}{{\int }_{-\infty }^{\infty }{e}^{-{\mathrm{Pe}}_{t}W(z,\theta )}dz}\right],$$

$$\dfrac{\rho_{0,1}(z,\theta)}{\bar{\rho}_{0}}\simeq\dfrac{1}{2}\left[1-\mathrm{Pe}_{t}\Big{\{}W(z,\theta)+\dfrac{1}{2}\int_{-1}^{1}W(z,\theta)dz+\frac{2}{3}\frac{R}{d}\int d\theta\int_{-1}^{1}\omega(z,\theta)dz\Big{\}}\right]+\mathcal{O}(\mathrm{Pe}_{t}^{2})\,.$$

$$\rho (x,z,\theta )=p(x)\frac{g(x,z,\theta )}{{\int }_{-\infty }^{\infty }{\int }_{-\pi }^{\pi }g(x,z,\theta )dzd\theta }$$

$$g(x,z,\theta )=\frac{e^{-{\mathrm{Pe}}_{t}W(x,z,\theta )}}{{\int }_{-\infty }^{\infty }{e}^{-{\mathrm{Pe}}_{t}W(x,z,\theta )}dz}\exp \left[\frac{4}{3}\frac{R}{d}{\mathrm{Pe}}_t\int d\theta \frac{{\int }_{-\infty }^{\infty }\omega (x,z,\theta )e^{-{\mathrm{Pe}}_{t}W(x,z,\theta )}dz}{{\int }_{-\infty }^{\infty }{e}^{-{\mathrm{Pe}}_{t}W(x,z,\theta )}dz}\right]$$

$${\int }_{-\infty }^{\infty }{\int }_{-\pi }^{\pi }\left\{\frac{{\partial }^2}{\partial x^2}\rho +\frac{\partial }{\partial x}\left[\rho \frac{\partial }{\partial x}\beta C\right]-{\mathrm{Pe}}_t\frac{\partial }{\partial x}\left[\left(\sin \theta +v_x\right)\rho \right]\right\}dzd\theta =0,$$

$${\int }_{-\infty }^{\infty }{\int }_{-\pi }^{\pi }\left\{\frac{\partial }{\partial x}\rho +\rho \left[\frac{\partial }{\partial x}\beta C-{\mathrm{Pe}}_t\left[\left(\sin (\theta )+v_x\right)\right]\right]\right\}dzd\theta =-I_x,$$

$${\int }_{-\infty }^{\infty }{\int }_{-\pi }^{\pi }\frac{\partial }{\partial x}\rho (x,z,\theta )dzd\theta =\frac{\partial }{\partial x}p(x)$$

$$v_D(x)=-\beta \frac{\partial }{\partial x}\chi (x)={\int }_{-\infty }^{\infty }{\int }_{-\pi }^{\pi }\left[-\frac{\partial }{\partial x}\beta C+{\mathrm{Pe}}_t\left(\sin (\theta )+v_x\right)\right]\frac{g(x,z,\theta )}{{\int }_{-\infty }^{\infty }{\int }_{-\pi }^{\pi }g(x,z,\theta )dzd\theta }dzd\theta ,$$

$$\frac{\partial }{\partial x}p(x)+\beta p(x)\frac{\partial }{\partial x}\chi (x)=-I_x,$$

$$p(x)=e^{-\beta \chi (x)}\left[-I_x{\int }_0^x{e}^{\beta \chi (x^{\prime })}d{x}^{\prime }+\Theta \right].$$

$$p(x)=\Theta e^{-\beta \chi (x)}.$$

$$v_d=I_{x}L.$$

$$-\beta \frac{\partial }{\partial x}\chi (x)={\mathrm{Pe}}_t{\int }_{-\infty }^{\infty }{\int }_{-\pi }^{\pi }\frac{\left(\sin (\theta )+v_x\right)g(x,z,\theta )}{{\int }_{-\infty }^{\infty }{\int }_{-\pi }^{\pi }g(x,z,\theta )dzd\theta }dzd\theta -\beta \frac{\partial }{\partial x}\mathcal{A}(x).$$

Full text

I Introduction

Organisms as well as synthetic particles swimming at low Reynolds number attain net displacement by locally stirring the fluid Lauga and Powers (2009). Therefore, their swimming performances are affected by the presence of boundaries, interfaces, or other particles, which perturb the fluid flow generated by the microswimmers and thereby ultimately influence their swimming speed. For example, sperm cells as well as bacteria, have been shown to accumulate at solid walls Rotschild (1963); Di Luzio et al. (2005); Lauga et al. (2006); Friedrich et al. (2010). Moreover, in the presence of solid boundaries or fluid-fluid interfaces E.Coli bacteria swim along circular trajectories Berg and Turner (1990); Di Leonardo et al. (2011); Pimponi et al. (2016). In many biological situations as well as technological applications microswimmers are required to move in confined regions as it happens in microfluidic devices or in the female reproductive tract Suarez and Pacey (2006). Recently, such situations have been addressed more systematically both in theory Popescu et al. (2009); Spagnolie and Lauga (2012); Elgeti et al. (2015); Uspal et al. (2015); Colberg and Kapral (2015); Wu et al. (2015); Domínguez et al. (2016a, b); Zöttl and Stark (2016); Malgaretti et al. (2016a); Menzel et al. (2016) as well as experiments Hulme et al. (2008); Drescher et al. (2011); Denissenko et al. (2012); Altshuler et al. (2013).

Up to now, much attention has been payed to the case in which active particles are confined by walls or interfaces. However, for passive systems it is well known that the shape of the boundaries can induce novel dynamic regimes absent in the case of homogeneous confinement Hänggi and Marchesoni (2009); Burada et al. (2009); Malgaretti et al. (2013a). Similarly, recent works have shown that active particles moving in a non-uniform geometrical confinement give rise to novel dynamic effects. For example, in experiments the motion of bacteria was rectified by funnel Galajda et al. (2007) or ratchet-like Di Leonardo et al. (2010) potentials. Moreover, theoretical models dealing with active particles under the influence of period enthalpic Pototsky et al. (2013) or entropic Ghosh et al. (2013); Ghosh (2014); Wu et al. (2015) potentials, show rectification of active motion under suitable conditions. Or, they demonstrate how active Brownian and run-and-tumble particles separate in circular mazes Khatami et al. (2016). Finally, dynamic collective effects arise, for example, when bacteria are confined in circular pools, where they spontaneously circle in one direction Lushi et al. (2014), or when active particles are confined in narrow channels joining larger reservoirs, where self-sustained density oscillations evolve Paoluzzi et al. (2015).

In this contribution, we study the dynamics of microswimmers moving in a channel with varying cross section, where channel walls are corrugated. In particular, approximating the fluid flow, initiated by a microswimmer, by the leading far–field expression Drescher et al. (2011); Spagnolie and Lauga (2012); Schaar et al. (2015), we derive an effective interaction between microswimmers and the corrugated channel walls that depends on both the swimming mechanism and the wall geometry. In order to derive analytical expressions, we assume that the channel cross section varies smoothly enough such that swimmers explore the cross-sectional direction while moving along the channel axis. In such a regime, the Fick–Jacobs approximation Zwanzig (1992); Reguera and Rubi (2001); Berezhkovskii et al. (2007); Kalinay and Percus (2008); Berezhkovskii and Szabo (2011); Dagdug and Pineda (2012); Chacón-Acosta et al. (2013) properly captures the overall dynamics of systems as diverse as molecular motors Malgaretti et al. (2012, 2013b, 2014a), driven particles Reguera et al. (2006, 2012), electrokinetic phenomena Malgaretti et al. (2014b, 2015, 2016b), and polymer translocation Bianco and Malgaretti (2016) just to mention a few among others (see recent review articles Hänggi and Marchesoni (2009); Burada et al. (2009); Malgaretti et al. (2013a) for a more comprehensive list).

Here, we extend the Fick–Jacobs approximation to active particles111A first attempt to generalize the Fick–Jacobs equation to the case of active Brownian particles (i.e. in the absence of hydrodynamic coupling) can be found in Ref. Sandoval and Dagdug (2014) and we will characterize the motion of microswimmers in narrow channels with corrugated walls. In particular, our model shows that the geometrical confinement and the hydrodynamic interactions of the microswimmer with the channel walls can induce opposing effective forces on the swimmer. Hence, there is a cross-over value of the Péclet number, for which these two contributions are equal in magnitude. The whole formalism is valid independent of a special swimmer type. In order to quantify the effective coupling between the geometrical confinement and the microswimmers, we use generic expressions of the flow field induced by model microswimmers such as pushers, pullers, and source dipoles similar to the approaches in Refs. Hernandez-Ortiz et al. (2005); Trouilloud et al. (2008); Baskaran and Marchetti (2009); Spagnolie and Lauga (2012); Lopez and Lauga (2014); Hennes et al. (2014); Schaar et al. (2015). The no-slip boundary condition at the channel walls is accounted for by the method of images Spagnolie and Lauga (2012). Our results show that the strength of accumulation of the swimmers at the channel walls is strongly affected by both the channel geometry and the underlying swimming mechanism of pushers, pullers, and source dipoles. In particular, fast pullers accumulate at different locations along the channel than pushers and source dipoles, hence allowing for a mechanism to separate active particles. Finally, for asymmetric channels we show that the motion of the swimmers becomes rectified and that different swimmer types can move in opposite directions along the channel axis.

The article is organized as follows. In section II we derive general expressions for our model of confined microswimmers. In sections II.A and II.B we then take the limit of small ratio of rotational to translational Péclet number, which means narrow channels, and apply our model to the two cases of plane walls and smoothly corrugated walls. In order to exploit the predictions of the model, in section III we use the generic expressions for the flow fields of model microswimmers to discuss their dynamics in channels with plane and corrugated walls. Finally, in section IV we summarize our results.

II Model

We analyze the dynamics of spherical swimmers moving between two corrugated walls, whose distance $2h(x)$ varies solely along the $x$ direction and is constant along the $y$ direction (see Fig. 1). In order to keep our analysis simple and to be able to formulate analytical insight, we assume that the swimming axes of the active particles is constrained to the $x-z$ plane222Experimentally such a configuration can be attained by trapping the swimmers by means of optical tweezers. and that the system extends infinitely along the $y$ direction. Such an assumption together with the homogeneity of the channel cross section along the $y$ axis makes the probability distribution independent of $y$ and the problem becomes two-dimensional. Accordingly, in the overdamped regime, the dynamics of a dilute suspension of non-interacting microswimmers is captured by the Smoluchowski equation that governs the temporal evolution of the density $\rho(\mathbf{r},\theta)$ of swimmers located at $\mathbf{\mathbf{r}}=(x,z)$ and with their swimming axes oriented at an angle $\theta$ against the $z$ axis Elgeti and Gompper (2013):

[TABLE]

where

[TABLE]

with

[TABLE]

Here, $\mathbf{J}_{\mathbf{r}}$ and $J_{\theta}$ are the respective translational and rotational probability current densities, which consist of diffusional ($D_{t}$, $D_{r}$), active ($v_{0}\mathbf{n}$), drift ($-\beta D_{t}\vec{\nabla}_{\mathbf{r}}C$), and wall-induced ($v_{0}\mathbf{v}$, $\omega_{0}\omega$) contributions. In particular, $D_{t}$ ($D_{r}$) is the translational (rotational) diffusion coefficient333We remark that when a particle approaches a wall, its translational and rotational diffusion coefficients become anisotropic and have to be described by tensors, the coefficients of which depend on the distance to the bounding wall. The anisotropy is due to the different hydrodynamic coupling perpendicular and parallel to the wall. Furthermore, one also has a wall-induced rotation-translation coupling tensor. In the following we are interested in the regime of moderate Péclet numbers with $\mathrm{Pe}_{t}>1$, for which diffusion is subdominant as compared to the active displacement of the microswimmers. Accordingly, in order to keep the derivation of the model as simple as possible, we neglect these corrections in the diffusion tensors and we regard $D_{t}$ and $D_{r}$ as constant scalar parameters., $v_{0}\mathbf{n}$ is the swimming velocity directed along the unit vector $\mathbf{n}$, and $\beta^{-1}=k_{B}T$ with $k_{B}$ the Boltzmann constant and $T$ the absolute temperature. The dimensionless quantities $\omega$ and $\mathbf{v}$ are the respective angular and linear velocities induced by the hydrodynamic coupling of the swimmer flow field with the channel walls (see appendix B for explicit expressions) and $\omega_{0}=v_{0}/R$ is the relevant angular-velocity scale with particle radius $R$. We remark that for active Brownian particles, where any hydrodynamic interactions with the walls are neglected, one has $\omega=0$ and $\mathbf{v}=\mathbf{0}$, while for passive particles, for which $v_{0}=\omega_{0}=0$, Eqs. (3) and (4) reduce to the standard diffusion equations including a translational drift current induced by the geometrical confinement.

The force potential $C$ in Eq. (3),

[TABLE]

confines particles between the channel walls located at $\pm h(x)$444Eq. (5) has been obtained by expanding the exact expression for the available space a rigid sphere of radius $R$ can explore in a channel with varying cross section about a channel with constant cross section, $h(x)=h_{0}$, assuming small values of $\partial_{x}h(x)$..

The steady-state probability density is obtained by solving

[TABLE]

which gives

[TABLE]

where we normalized all lengths by $d=h_{0}-R$. Furthermore, we identified the rotational, $\mathrm{Pe}_{r}=\frac{v_{0}}{D_{r}d}$, and translational, $\mathrm{Pe}_{t}=\frac{v_{0}d}{D_{t}}$, Péclet numbers. Equation (7) is complemented by the following boundary conditions using the unit vector $\mathbf{n}_{h}$ for the normal at the channel walls:

[TABLE]

where $L$ is the length of the channel. We remark that Eq. (12) stems from the mirror symmetry about the channel axis, which implies $\rho(x,z,\theta)=\rho(x,-z,\pi-\theta)$ and $\omega(x,z,\theta)=-\omega(x,-z,\pi-\theta)$. This implies that the integrated rotational flux as formulated in Eq. (12) vanishes for all positions along the channel axis hence preventing the onset of local recirculation of microswimmers.

Equation (7) is quite involved because it accounts for contributions stemming from both the translational ($\mathbf{J}_{\mathbf{r}}$) and rotational ($J_{\theta}$) fluxes. In order to disentangle these contributions, we will follow Ref. Elgeti and Gompper (2013) and assume $\mathrm{Pe}_{t}/\mathrm{Pe}_{r}\ll 1$. We remark that for the thermal or Stokes-Einstein values of the diffusion coefficients $D_{t}$ and $D_{r}$ the condition $\mathrm{Pe}_{t}/\mathrm{Pe}_{r}=3d^{2}/4R^{2}\ll 1$ gives $d^{2}\ll R^{2}$, which implies narrow channels. The ratio $\mathrm{Pe}_{t}/\mathrm{Pe}_{r}$ can be affected by the activity of the microswimmers. For example, for a self-diffusiophoretic particle with $R=1\mu\text{m}$ the rotational decorrelation time has been measured Howse et al. (2007) to be $D^{-1}_{r}\sim 4\,\text{sec}$. Hence, the condition $\mathrm{Pe}_{t}/\mathrm{Pe}_{r}=d^{2}D_{r}/D_{t}\ll 1$ together with $D_{t}=k_{B}T/(6\pi\eta R)$ gives an upper bound $d\simeq 1\mu\text{m}$, i.e., the channel half width is twice the particle radius. Finally, we note that the regime $d^{2}\ll R^{2}$ can be attained for all synthetic and biological swimmers by properly tuning the accessible average width of the channel, $2d$. We now formally solve Eq. (7) for the two cases of plane and corrugated channel walls using an expansion in $\mathrm{Pe}_{t}/\mathrm{Pe}_{r}\ll 1$.

II.1 Plane channel walls

For a channel with plane walls, $h(x)=h_{0}$, the probability density $\rho$ and the velocities $v$, $\omega$ are uniform along $x$. In the absence of any spontaneous symmetry breaking along the $x$ axis, this gives $J_{x}=0$. The boundary condition, Eq. (8), reads:

[TABLE]

In zeroth order in $\mathrm{Pe}_{t}/\mathrm{Pe}_{r}$ the rotational contribution in Eq. (7) is negligible and we have:

[TABLE]

Since $\rho$ does not depend on $x$, which means zero flow along $x$ ($J_{x}=0$), Eq. (14) can be integrated once and reduces to

[TABLE]

The constant probability flux $J_{z,0}(\theta)$ vanishes due to Eq. (13). Equation (15) can be integrated and one obtains

[TABLE]

In the last equation we have introduced the effective potential

[TABLE]

using the properties of the hard-core potential $C(z)$ from Eq. (5). We recall that all lengths are measured in units of $d=h-R$. The factor $\lambda_{0}$ is determined by requiring that $J_{\theta}$, once calculated using Eq. (16), fulfills the local isotropy of Eq. (12). In particular, if $J_{x}=0$, Eq. (12) reduces to555 Equation (18) is obtained by integrating the steady-state condition along $z$, namely $\int dz(\partial_{x}J_{x}+\partial_{z}J_{z}+\partial_{\theta}J_{\theta})=0$. Then, using Eq. (13) and the fact that for plane channel walls we have $J_{x}=0$, this condition reduces to $\int dz\partial_{\theta}J_{\theta}=\partial_{\theta}\int dzJ_{\theta}=0$, which implies that $\int dzJ_{\theta}$ is a constant. Finally, due to Eq. (12) we obtain $\int dzJ_{\theta}=0$.

[TABLE]

Rewriting $\rho_{0}=\lambda_{0}\hat{\rho}_{0}$, using $J_{\theta}$ from Eq. (4), and Eq. (18) we obtain:

[TABLE]

whose solution reads:

[TABLE]

Finally, we substitute the last expression in Eq. (16) and write

[TABLE]

where $\bar{\rho}_{0}$ is determined by the normalization condition $\int_{-\pi}^{\pi}d\theta\int_{-\infty}^{\infty}\rho_{0}(z,\theta)dz=1$. We also rewrote $\omega_{0}/D_{r}=4/3v_{0}R/D_{t}=4/3\mathrm{Pe}_{t}R/d$ using the thermal values for $D_{r}$ and $D_{t}$. We remark that neglecting the hydrodynamic coupling with the walls, namely for $v=\omega=0$, Eq. (21) reduces to Eq. (9) of Ref. Elgeti and Gompper (2013). Interestingly, the probability distribution of Eq. (21) does not factorize in the two “natural” effective potentials, namely $W(z,\theta)$ and $\int\omega(z,\theta)d\theta$, from which the stall force and torque, necessary to stop the microswimmer, can be derived. This implies that the translational and rotational degrees of freedom are quite entangled. In order to get insight into their relative contributions to the particle probability distribution, we expand Eq. (21) for small values of $\mathrm{Pe}_{t}$ and for $|z|<1$. In particular, at first order in $\mathrm{Pe}_{t}$ we obtain:

[TABLE]

Note that the angular velocity only contributes with its mean value and this term vanishes for pure source-dipole swimmers [see Eq. (63), which shows that $\omega$ is an odd function in $z$.]

II.2 Corrugated channel walls

In order to characterize the dynamics along the channel axis, we determine and discuss the probability density $p(x)$ to find the microswimmer at position $x$ for any $\theta$, $z$ and also the mean drift velocity in $x$ direction. To do so, we take inspiration from the Fick–Jacobs approximation and extend it to the case of active particles. Accordingly, we start with an ansatz for the full probability density:

[TABLE]

where

[TABLE]

has the same functional form as in Eq. (21) for plane channel walls and thus fulfills the boundary conditions of Eqs. (11), (12), and (8) (for the latter see Appendix A of Ref.Marini Bettolo Marconi et al. (2015)). In making this ansatz we are again in the limit of $\mathrm{Pe}_{t}/\mathrm{Pe}_{r}\ll 1$ and also assume small variations of the channel width, $|\partial h(x)/\partial x|\ll 1$. We also note that $\omega$ and $v$ account for the local slope of the channel walls and thus they now depend on $x$ (see Appendix B). This ansatz for $\rho(x,z,\theta)$ in Eqs. (23) and (24) represents the extension of the Fick–Jacobs approximation to active particles.

We integrate Eq. (7) along $z$ and $\theta$ and use the boundary conditions from Eqs. (8),(12) together with $|\partial_{x}h(x)|\ll 1$ Marini Bettolo Marconi et al. (2015) to obtain

[TABLE]

or, equivalently, after integrating once:

[TABLE]

where $I_{x}$ is the probability current along the $x$ direction. We use the ansatz from Eq. (23) to show that

[TABLE]

and also define the average drift velocity $v_{D}(x)$ at axial position $x$,

[TABLE]

as the derivative of the effective total potential $\chi(x)$, which in one dimension can always be defined. Using Eqs. (27) and (28), we are able to write Eq. (26) in compact form,

[TABLE]

the solution of which reads

[TABLE]

Here, $\Theta$ and $I_{x}$ are determined by imposing the boundary condition of Eq. (9), and the normalization $\int p(x)dx=1$. For $I_{x}=0$, Eq. (30) reduces to

[TABLE]

Finally, we can use the current $I_{x}$ to define the net drift velocity along the $x$ axis as

[TABLE]

This completes the full solution of our problem. With the help of $\rho(x,z,\theta)$ from Eqs. (23) and (24), we can calculate $v_{D}(x)$ from Eq. (28), integrate it once to obtain $\chi(x)$, and then determine the probability density $p(x)$ from Eq. (30).

Even though the functional form of Eq. (30) is similar to that obtained for passive systems governed solely by conservative forces Reguera and Rubi (2001), the difference between Eq. (30) and the corresponding expression in the passive case lies in the form of the effective total potential, $\chi(x)$, the derivative of which can be expressed as

[TABLE]

Here, we have introduced the gradient of an effective entropic potential $\mathcal{A}(x)$, which generalizes the entropic potential of passive systems Reguera and Rubi (2001) and derives from

[TABLE]

To arrive at the last expression, we used the explicit form of $C(x,z)$ from Eq. (5) (see Appendix A for details). In particular, for passive systems with only conservative interactions with bounding walls, $\mathrm{Pe}_{t}=0$ and from Eq. (26) one has $g(x,z,\theta)\propto\exp[-\beta C(x,z)]$. Accordingly, the first term on the right-hand side of Eq. (33) vanishes and Eq. (34) can be integrated leading to

[TABLE]

Therefore, for passive systems, $\mathcal{A}(x)$ reduces to the usual entropic potential666The term “entropic” stems from the fact that $h(x)-R$ measures the number of states available in the transverse direction and therefore $\ln[(h(x)-R)/d]$ is proportional to the local entropy. of Eq. (35). In contrast, in the present active case $\mathrm{Pe}_{t}\neq 0$ and the effective interactions are not conservative. Accordingly, $g(x,z,\theta)\nsim\exp[-\beta C(x,z)]$, and the form of both $\chi(x)$ and $\mathcal{A}(x)$ is more involved. Interestingly, Eq. (33) identifies two separate contributions to the mean drift velocity $v_{D}(x)$, namely an “active” one [the first term on the rhs of Eq. (33)] and an “entropic” one [the second term on the rhs of Eq. (33)]. Therefore, by tuning the geometry of the channel and the activity of the particles via the Péclet number, it is possible to tune the relevance of the entropic against active contribution. Moreover, the sign of the active contribution depends on the hydrodynamic coupling to the bounding walls and hence on the characteristic flow profile induced by the swimmer. Accordingly, for those cases for which the entropic and the active contribution have opposite sign, Eq. (33) predicts the existence of a critical value of the Péclet number, $\mathrm{Pe}_{c}$, for which the dynamics crosses over from entropy-controlled (small $\mathrm{Pe}_{t}$) to activity-controlled (large $\mathrm{Pe}_{t}$). In order to gain some insight into the magnitude of $\mathrm{Pe}_{c}$, we assume

[TABLE]

namely the density is constant inside the channel and zero outside777The probability distribution in Eq. (36) is the one that should be expected for passive systems, i.e., for $\mathrm{Pe}_{t}=0$. Here, we assume that Eq. (36) roughly captures the functional form of $g(x,z,\theta)$ also for $\mathrm{Pe}_{t}\neq 0$. As discussed in the Results section, such an assumption is justified a posteriori due to the good matching of the predictions that can be derived, Eq. (39), when compared with the numerical integration of $g(x,z,\theta)$.. In this case, using Eq. (5), Eq. (33) reads:

[TABLE]

where $\left<v_{x}\right>=\int_{-\infty}^{\infty}\int_{-\pi}^{\pi}v_{x}\bar{g}dzd\theta$ and $\int_{-\infty}^{\infty}\int_{-\pi}^{\pi}\sin(\theta)\bar{g}dzd\theta=0$ due to the symmetry of $\bar{g}$. Setting $v_{D}(x)=0$ in Eq. (37), we identify the cross-over value of the Péclet number,

[TABLE]

where in the last step we used the amplitude $h_{1}$ of the corrugation of the channel walls and we made the approximation $\partial_{x}h(x)\simeq\frac{2h_{1}}{L/2}$. Moreover, assuming that the magnitude of $\left\langle v_{x}\right\rangle$ is comparable to $v_{0}$, i.e., for $\left\langle v_{x}\right\rangle\simeq 1$ we obtain

[TABLE]

III Results

As shown by Eqs. (21) and (30), the dynamics of an active particle inside a channel is strongly affected by the hydrodynamic coupling to the bounding walls. Since an exact solution of the Stokes equation for such a geometry is lacking, we have to rely on approximate solutions or generic models. In particular, in order to derive general results that do not depend on a specific swimming mechanism, we study model microswimmers such as pushers, pullers, and source dipoles, the velocity fields of which correspond to the lowest-order contributions in a far–field multipole expansion of real microswimmers. Moreover, we neglect that their finite size disturbs the fluid velocity field of the generic swimmers as it was done in other studies Hernandez-Ortiz et al. (2005); Trouilloud et al. (2008); Baskaran and Marchetti (2009); Spagnolie and Lauga (2012); Lopez and Lauga (2014); Hennes et al. (2014); Schaar et al. (2015). For such model microswimmers the effective interactions with the channel walls can be captured by the method of images888In order to keep our treatment simple, we have truncated the set of images at the first iteration, see Appendix B.. Interestingly, it has been shown that for some cases Spagnolie and Lauga (2012) these generic far–field expansions quantitatively capture the hydrodynamic interactions with bounding walls even when the microswimmers are very close to these walls. As a result, hydrodynamic interactions between model microswimmers and channel walls are accounted for by the linear and angular velocities induced by the image flow fields at the center of mass of the microswimmers. Finally, the fact that real microswimmers have a finite extent and in order to avoid divergences in the hydrodynamic coupling, we prevent model microswimmers from approaching channel walls not closer than their “radius” $R$.

In this perspective, we exploit generic far–field approximations to grasp the relevance of the swimmer type and mechanism, and we characterize the effective dynamics of four different kinds of swimmers Spagnolie and Lauga (2012). In particular, we consider swimmers, which generate either a source dipole of strength $q$ or a force dipole of strength $p$, distinguishing pushers ($p>0$) and pullers ($p<0$), and compare them to pure active Brownian particles, which do not hydrodynamically couple to the bounding walls. The different swimmer types are determined by the wall-induced linear ($\mathbf{v}$) and angular ($\omega$) velocity fields, for which we take the far–field expressions as summarized in Appendix B. Source dipoles also called “neutral swimmers”, the flow fields of which decay as $1/r^{3}$, resemble the motion of biological microswimmers such as Paramecia Ishikawa and Hota (2006) or artificial swimmers as the active emulsion droplets of Refs. Thutupalli et al. (2011); Schmitt and Stark (2013). Typical force dipoles in nature are bacteria such as E.Coli that are characterized by $p>0$ and are named “pushers” since they push the fluid along the direction of motion, whereas algae like Clamydomonas are characterized by $p<0$ and are named “pullers” since they pull on the fluid along the direction of motion. Their flow fields always decay as $1/r^{2}$. In the following we characterize the steady-state probability distribution along a microchannel with periodic boundary conditions at both channel ends.

Plane channel walls

Substituting the explicit expression for the angular [Eq. (63)] and linear [Eq. (64)] velocity provided in Appendix B into Eq. (21), we can calculate the density profile in steady state. At first, we focus on the lateral probability density,

[TABLE]

which for plane walls does not depend on the axial coordinate $x$. As shown in Fig. 2(a) (and in agreement with Ref. Elgeti and Gompper (2013)) an active Brownian particle shows an excess accumulation at the channel walls as compared to a passive particle embedded in the same channel. When an active particle stirs the fluid, the induced hydrodynamic coupling alters the previous picture, as shown in Fig. 2(a). In particular, while pullers [orange long-dashed line in Fig. 2(a)] show an enhanced accumulation at channel walls as compared to active Brownian particles, pushers [blue dashed line in Fig. 2.(a)] and neutral swimmers [red dot-dashed line in Fig. 2.(a)] experience a decreased accumulation at the channel walls. This is in agreement with the larger detention times of pullers found in Ref. Schaar et al. (2015) compared to active Brownian particles. Interestingly, the accumulation of active particles at channel walls, $\rho_{\theta}(z=1)$, shows a non trivial dependence on $\mathrm{Pe}_{t}$, as shown in Fig.2.(b). Looking at the dependence of $\rho_{\theta}(z=1)$ upon increasing $\mathrm{Pe}_{t}$, we see that active Brownian particles, pushers, and pullers are characterized by a monotonic increase of $\rho_{\theta}(z=1)$. In contrast, the density of neutral swimmers the walls shows a plateau for Péclet values larger than $\mathrm{Pe}_{t}\simeq 10$.

The diverse accumulation at channel walls is related to the different angular orientations at contact. In fact, as shown in Fig. 2.(c), active Brownian particles and pullers accumulate at the walls mainly orthogonally pointing towards the wall. In contrast, pushers and neutral swimmers tend to align parallel with the walls (see the peaks of $\rho(z=\pm 1,\theta)$ close to $\theta=\pi/2$). Finally, while the orientation of active Brownian particles and pullers does not depend on the value of $\mathrm{Pe}_{t}$, pushers and neutral swimmers show a non trivial dependence on $\mathrm{Pe}_{t}$ and they tend to fully align parallel to the channel walls for increasing values of $\mathrm{Pe}_{t}$, as shown in Fig. 2.(d). On the contrary, for $\mathrm{Pe}_{t}\lesssim 1$, the angle $\theta_{\text{max}}$, where $\rho(z=1,\theta)$ attains its maximum, is captured by Eq. (22). For the case of a plane channel it reads

[TABLE]

where we have used Eqs. (63),(64). In particular, for active Brownian particles ($q=p=0$) Eq. (41) reduces to

[TABLE]

whereas for neutral swimmers ($q\neq 0$, $p=0$) and pushers/pullers ($q=0$, $p\neq 0$) Eq. (41) becomes

[TABLE]

Interestingly, while active Brownian particles always point mainly towards the wall [$\theta_{\text{max}}=0$ for Eq. (42)], the maximum of $\rho_{0,1}(z=1,\theta)$ for neutral and force-dipole swimmers [see Eqs. (43) and (44), respectively] has a more involved parameter dependence. In particular, for $q\frac{h_{0}}{(1-h_{0}^{2})^{2}}<\frac{1}{2}$ [valid for the values of $h_{0}$, $d$, and $q$ used in Fig. 2 (d)] the density in Eq. (43) is maximized for $\theta_{\text{max}}=0$, whereas for $q\frac{h_{0}}{(1-h_{0}^{2})^{2}}>\frac{1}{2}$ one has $\theta_{\text{max}}=\pi$. For pushers/pullers the density in Eq. (44) has quite an involved dependence on the parameters. In particular, for the geometric values used in Fig. 2, it has its maximum at $\theta_{\text{max}}=0$ for pullers ($p<0$), whereas for pushers ($p>0$) one obtains $\theta_{\text{max}}\simeq\pi/5$, which is the value, $\theta_{\mathrm{max}}$ tends to in Fig.  2(d).

Corrugated channel walls

We have more thoroughly characterized the dynamics of active particles in a channel with corrugated walls, the half width of which is defined as

[TABLE]

Here, $2h_{0}$ is the average channel width, $h_{1}$ the amplitude of the wall modulation, and $L_{1}$ characterizes the asymmetry of the wall corrugation along the channel axis (see Fig. 1). In the case of modulated channel walls the explicit expression of the linear ($\mathbf{v}$) and angular ($\omega$) velocity, which we use explicitly in Eqs. (23) and (24) and all the following implications, must take into account the modulation of the channel width, $2h(x)$. If the channel cross section varies on lengths much larger than particle size, $L\gg R$, locally we can regard channel walls as infinite planes tilted by an angle $\psi(x)=\arctan(\partial_{x}h(x))$ with respect to the channel axis. For this case, the hydrodynamic coupling between our model microswimmers and the tilted channel walls is known in the far–field regime Spagnolie and Lauga (2012) (see Appendix B). For corrugated channel walls, the mean drift velocity $v_{D}(x)$ along the channel axis [see Eq. (33)] is generally non-vanishing and the local drift might result in a net rectification of particle flow giving a net current $I_{x}$ or a net drift velocity $v_{d}=I_{x}L$, as defined in Eq. (32). Substituting the expressions for the linear velocity along the normal direction [$v_{z}(x,z,\theta)$ from Eq. (61)], along the longitudinal or axial direction [$v_{x}(x,z,\theta)$ from Eq. (62)], and the angular velocity [$\omega(x,z,\theta)$ from Eq. (60)] into Eqs. (17), (24), and (28), we can calculate the probability distribution along the axial position from Eq. (30).

For $L_{1}=L/2$ the channel corrugation is symmetric hence the net current is zero, $I_{x}=0$. Figure 3 shows the behavior of $p(x)$, as obtained from Eq. (31), as a function of the axial position $x$. Interestingly, Fig. 3(a) already shows a significant deviation in the axial density of an active Brownian particle compared to a passive one. All particles accumulate more at positions, where the channel is broadest, and the excess of accumulation increases strongly for larger values of $\mathrm{Pe}_{t}$. Looking at panels (b), (c) and (d) of Fig. 3, we recognize that hydrodynamic coupling strongly modulates the distribution of swimmers as compared to active Brownian particles. In particular, neutral swimmers [panel (b)] and pushers [panel (c)] behave qualitatively similar to active Brownian particles, but the net effect of the hydrodynamic coupling induced by the walls is to enhance particle accumulation, where the channel is broadest. On the contrary, pullers [panel (d)] have quite a different behavior. Panel (d) of Fig. 3 shows that for smaller values of $\mathrm{Pe}_{t}$ pullers like the previous swimmers accumulate more in the broader parts of the channel. However, for increasing values of $\mathrm{Pe}_{t}$ the puller density becomes more uniform ($\mathrm{Pe}_{t}\simeq 3$) and, eventually, for $\mathrm{Pe}_{t}\geq 5$, very interestingly, pullers accumulate at channel bottlenecks.

Figure 4(a) shows in detail the maximum of the probability distribution, $p_{max}$, as a function of $\mathrm{Pe}_{t}$. Clearly, for active Brownian particles, neutral swimmers, and pushers, the maximum $p_{max}$ increases monotonically upon increasing $\mathrm{Pe}_{t}$ and is always located at $x=L/2$. On the contrary, pullers are characterized by a “critical” value of the Péclet number, $\mathrm{Pe}_{c}$, below which the maximum is located at $x=L/2$, whereas for $\mathrm{Pe}_{t}>\mathrm{Pe}_{c}$ the maximum is located at $x=0,L$. Interestingly, the value of $\mathrm{Pe}_{c}$ is properly captured by Eq. (39), which indicates the Péclet number, where the swimmer dynamics crosses over from entropy-controlled to activity-controlled. It gives $\mathrm{Pe}_{c}\simeq 2$ for $h_{1}=0.4$ and $L=10$, which are the parameters used to obtain the data in Figs. 3 and 4. Finally, we note that for pushers and neutral swimmers the value of $p_{max}$ is quite independent of the channel corrugation, $h_{1}$, whereas pullers and active Brownian particles are sensitive to $h_{1}$, as Fig. 4(b) indicates.

When the channel is not symmetric, i.e., $L_{1}\neq L/2$ in Eq. (47), in principle a non-zero net drift velocity $v_{d}$ [see Eq. (33)] occurs meaning rectification of the swimmer flow. Therefore, the probability distribution $p(x)$ has to be calculated from the more general expression of Eq. (30). In Fig. 5 we plot it for a range of $\mathrm{Pe}_{t}$ values for different swimmer types. The different graphs show that active Brownian particles as well as the other swimmer types behave similar as in a symmetric channel concerning the location of maximum accumulation and the variation with $\mathrm{Pe}_{t}$. However, the probability distribution $p(x)$ reflects the asymmetry of the channel and gives a non-zero net drift. Indeed, as demonstrated in Fig. 6, the different hydrodynamic swimmer types show a non-zero drift velocity $v_{d}$, while for active Brownian particles, which do not interact hydrodynamically with the channel walls, $v_{d}=0$. In particular, we find that both the sign and the magnitude of $v_{d}$ depends on the underlying physical mechanism leading to active displacement (hydrodynamic swimmer type) as well as on the value of $\mathrm{Pe}_{t}$ in a non-trivial way. Notably, we observe a non-monotonous behavior for all swimmers and for pullers also an inversion in the velocity. Finally the net swimmer flux vanishes if the channel becomes symmetric ($L_{1}=L/2$) as demonstrated in Fig. 6(b).

IV Conclusions

We have studied the dynamics of microswimmers in channels with corrugated walls. Under the assumption that the cross section of the channel varies smoothly, we have extended the Fick–Jacobs approximation to the case of active particles. Accordingly, the activity and the hydrodynamic coupling between the microswimmers and the channel walls can be captured by an effective potential $\chi(x)$ defined in Eq. (33). Interestingly, our model identifies two contributions to $\chi$, an active one proportional to $\mathrm{Pe}_{t}$ [first term on the right-hand side of Eq. (33)], and an entropic one [second term in Eq. (33)]. Unlike passive systems, the effective interaction between the microswimmers and the channel walls cannot be formulated with a conservative potential and, therefore, the entropic contribution to $\chi$ cannot be expressed in the form of Eq. (35) familiar from passive systems. In particular, our framework identifies a cross-over value, $\mathrm{Pe}_{c}$, of the Péclet number that separates two distinct regimes. For $\mathrm{Pe}_{t}\ll\mathrm{Pe}_{c}$ the system is in an entropy–controlled regime, in which the probability distribution is mainly determined by the geometry of the channel and thereby resembles the behavior of passive particles. On the contrary, for $\mathrm{Pe}_{t}\gg\mathrm{Pe}_{c}$ the system switches to an activity–controlled regime, in which the probability distribution depends strongly on the details of the swimming mechanism.

In order to quantify our findings, we have applied our framework to different types of model microswimmers by taking into account their hydrodynamic interactions with the channel walls using the far–field approximation of Refs. Happel and Brenner (1965); Spagnolie and Lauga (2012). However, we also disregarded that the finite size of microswimmers disturbs these generic velocity fields. Under these assumptions, we have been able to characterize the behavior of our model microswimmers in microchannels and to contrast it against active Brownian particles, i.e., active particles, which do not interact hydrodynamically with the walls. Our results show that the effective dynamics of confined microswimmers is strongly affected by their swimming mechanism. In fact, while pullers accumulate stronger than active Brownian particles at channel walls, for pushers and neutral swimmers accumulation is weaker. Such a behavior depending on the swimmer type can already be observed in channels with plane walls (see Fig. 2).

For corrugated channels the behavior of different generic microswimmer types is further amplified. Microswimmers show a non-uniform accumulation at the channel walls. Active Brownian particles, pushers, and neutral swimmers accumulate preferentially at positions, where the channel cross section is widest, irrespective of the value of $\mathrm{Pe}_{t}$. The accumulation is strongest for neutral swimmers followed by pushers and active Brownian particles. Pullers have a more involved dependence on $\mathrm{Pe}_{t}$. While for small values of $\mathrm{Pe}_{t}$ they also accumulate in the regions of widest cross section, for $\mathrm{Pe}_{t}>\mathrm{Pe}_{c}$ they preferentially accumulate at channel bottlenecks (see Fig. 3, Fig. 4), which is a clear signature of active motion. Finally, we have investigated the case, in which the channel corrugation breaks the fore-aft symmetry. This initiates a net flux of microswimmers along the channel axis, the direction of which strongly depends on both the underlying swimming mechanism and the value of $\mathrm{Pe}_{t}$ that controls the hydrodynamic interactions with the walls.

Our work clearly demonstrates how the hydrodynamics of different swimmer types determines their dynamics in channels with plane and corrugated channel walls. This has a clear potential for separating passive from active particles or different microswimmers from each other, which, for example, is important in biomedical applications.

Acknowledgments

P.M. acknowledges I. Pagonabarraga and P. Nowakowski for useful discussions and the Research Training Group GRK1558 funded by the Deutsche Forschungsgemeinschaft for travel grant and stay during scientific visiting at Technische Universität Berlin.

Appendix A Derivation of Eq. (34) of the main text

Here we derive the last term in Eq. (34). We divide the integral over $z$ into two parts:

[TABLE]

In the following we show how to rewrite the second term in the last expression. The first one can be treated following the same approach. In order to compute the derivative of the potential $C(x,z)$, it is useful to introduce the generalized Heaviside step function $T(x)$:

[TABLE]

where $a,b>0$ are small numbers999The choice of two different parameters, namely $a$ and $b$ for the regularization allows us to show that, as should be expected, the results do not depend on the value of the regularized distribution at $z=h(x)-R$., i.e. definitely smaller than $h(x)$. We remark that in the limit $a,b\rightarrow 0$ $T(x)$ recovers the Heaviside step function. Accordingly, we rewrite $C(x,y)$ in its regularized form:

[TABLE]

From the definition of the potential $W(x,z,\theta)$ in Eqs. (16),(17), we have:

[TABLE]

where $W_{0}(x,z,\theta)=-z\cos\theta-\int v_{z}(z,\theta)dz$. We remark that in the limit $\bar{c}\rightarrow\infty$ and $a,b\rightarrow 0$, Eqs. (50),(51) recover Eqs. (5),(17). According to Eq. (50) we have:

[TABLE]

Thus, the last integral in Eq. (48) reads:

[TABLE]

Before calculating the integral in the r.h.s of Eq. (53) we remark that in the limit $\bar{c}\rightarrow\infty$ the function $g(x,z,\theta)$ depends exponentially on $\bar{c}$ and therefore it is not possible to approximate the integral by the value of $g(x,z,\theta)$ at $z=h(x)-R$. Therefore, in order to calculate the integral in the r.h.s of Eq. (53) we isolate the contributions in $g(x,z,\theta)$ that depend on $\bar{c}$ and integrate them separately from those contribution that are not dependent on $\bar{c}$. Using Eq.(24) we rewrite $g(x,z,\theta)$ as

[TABLE]

Accordingly, using Eq. (51) the integral in Eq. (53) gives:

[TABLE]

For small values of $a$ and $b$ we have $e^{-\mathrm{Pe}_{t}W_{0}(x,z,\theta)}\simeq e^{-\mathrm{Pe}_{t}W_{0}(x,h(x)-R,\theta)}$. Hence, taking $e^{-\mathrm{Pe}_{t}W_{0}(x,h(x)-R,\theta)}$ out of the integral in $z$ and performing the integration the last expression can be approximated by

[TABLE]

where the last integral leads to:

[TABLE]

Finally substituting Eqs. (56),(57) in Eq. (48) we obtain:

[TABLE]

where we have used g(x,h(x)-R,\theta)=\tilde{g}(x,\theta)e^{-\mathrm{Pe}_{t}W_{0}(x,h(x)-R,{\color[rgb]{0,0,0}\theta})} and we have multiplied by $2$ since in Eq. (48) there are two contributions on the rhs that, due to the axial symmetry of the channel provide an equal contribution to Eq. (48). Finally, in the limit $\bar{c}\rightarrow\infty$ we obtain:

[TABLE]

Appendix B Derivation of the linear and angular velocity for channels with varying cross section

In order to account for the constraints imposed on the velocity field by the channel walls we exploit the method of images Happel and Brenner (1965); Spagnolie and Lauga (2012). In the presence of two walls in principle we should consider an infinite series of images since the we have to take into account that the images introduced to fulfill the no-slip boundary conditions on a wall will affect the value of the velocity field on the opposite wall. However, for the case under study at leading order in the far–field expansion the velocity field induced by microswimmers decays as $1/r^{2}$ for pushers/pullers and $1/r^{3}$ for neutral swimmers, where $r$ is the distance from the particle. In order to derive analytical expressions we disregard the contributions of higher-order images and we approximate the series with its first-order contribution.

Under these assumptions, in the case of a channel whose cross section varies on length scales much larger than particle size, we can locally approximate channel walls as flat planes tilted by an angle $\psi(x)=\arctan(\partial_{x}h(x))$ with respect to the channel longitudinal axis. According to Ref. Spagnolie and Lauga (2012), in dimensionless units, we have:

[TABLE]

where $p=p^{\prime}/d^{2}$ and $q=q^{\prime}/d^{3}$ are the dimensionless amplitudes of the force dipole and source dipole respectively. When the channel cross section is not constant, the normal at the channel walls are not aligned with the $z$ direction. Therefore, the linear velocities read Spagnolie and Lauga (2012):

[TABLE]

where $v_{z,1,2}$ and $v_{x,1,2}$ are the corrections to the transverse and longitudinal velocities due to the two walls. Furthermore, $\mathbf{e}_{z,1,2}$ and $\mathbf{e}_{x,1,2}$ are the unit vectors along the $z$ and $x$ axis in the frame of reference aligned along the channel walls (see Fig. 7), while $\mathbf{e}_{x,z}$ are the respective unit vectors along and perpendicular to the channel walls. In particular, we recall that

[TABLE]

Therefore we have:

[TABLE]

The magnitudes of the contributions are:

[TABLE]

Therefore the linear velocities read:

[TABLE]

and

[TABLE]

For plane channel walls we have $h(x)=h_{0}$ and $\psi(x)=0$. Accordingly Eqs. (60),(61),(62) reduce to:

[TABLE]

Appendix C Enhanced hydrodynamic interactions for plane channel walls

Eqs. (63)-(65) are derived by adding up the two contributions stemming from the hydrodynamic interactions with both walls. In order to discuss the relevance of these wall-induced velocities, we compare Eqs. (63)-(65) with the linear and angular velocity of a particle close to a single wall:

[TABLE]

In particular, we are interested in the ratio between the contribution to the linear and angular velocity due to the hydrodynamic interactions with one or two walls, namely:

[TABLE]

Interestingly, according to Eqs. (69)-(72) the values of $\delta_{p,q}$ and $\gamma_{p,q}$ do not depend on the orientation of the microswimmers for both pushers/pullers and neutral swimmers. As shown in Fig. 8, the two-walls scenario is quite different from the single-wall case. When particles are close to the wall (i.e. $z\rightarrow-1$ for the parameters of Fig. 8), the contribution of the second wall vanishes and the value of $\gamma_{p,q}$ and $\delta_{p,q}$ approach unity. In contrast, when particles are closer to the center of the channel (i.e. $z\rightarrow 0$) the presence of a second wall strongly affects the overall modulation in the linear and angular velocity as compared to the single wall case.