Controlling effective dispersion within a channel with flow and active walls

TL;DR

This paper explores how flow profiles and active or pulsating channel walls influence Taylor dispersion, providing analytic expressions to control solute transport in biological and artificial systems.

Contribution

It derives closed-form expressions for effective dispersion considering active walls and flow profiles, revealing regimes of dispersion enhancement or suppression.

Findings

Active walls can decrease or increase dispersion depending on conditions.

Dispersion regimes include entropic slowdown, shuttle dispersion, and Taylor dispersion.

Analytic formulas enable targeted control of solute transport.

Abstract

Channels are fundamental building blocks from biophysics to soft robotics, often used to transport or separate solutes. As solute particles inevitably transverse between streamlines along the channel by molecular diffusion, the effective diffusion of the solute along the channel is enhanced - an effect known as Taylor dispersion. Here, we investigate how the Taylor dispersion effect can be suppressed or enhanced in different settings. Specifically, we study the impact of flow profile and active or pulsating channel walls on Taylor dispersion. We derive closed analytic expressions for the effective dispersion equation in all considered scenarios providing hands-on effective dispersion parameters for a multitude of applications. In particular, we find that active channel walls may lead to three regimes of dispersion: either dispersion decrease by entropic slow down at small Peclet number, or dispersion increase at large Peclet number dominated either by shuttle dispersion or by Taylor dispersion. This improves our understanding of solute transport e.g. in biological active systems such as blood flow and opens a number of possibilities to control solute transport in artificial systems such as soft robotics.

Equations

$$\frac{\partial C}{\partial t}=-U\frac{\partial C}{\partial z}+\left(\kappa +\frac{U^2{a}^2}{48\kappa }\right)\frac{{\partial }^{2}C}{\partial z^2},$$

$$\frac{\partial c}{\partial t}=-u\frac{\partial c}{\partial z}-v\frac{\partial c}{\partial r}+\kappa \frac{{\partial }^{2}c}{\partial z^2}+\kappa \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial c}{\partial r}\right),$$

$$\left(\frac{L}{U}\frac{\kappa }{a^2}\right)\frac{1}{r^{\prime }}\frac{\partial }{\partial r^{\prime }}\left(r^{\prime }\frac{\partial c}{\partial r^{\prime }}\right)=\frac{\partial c}{\partial t^{\prime }}+u^{\prime }\frac{\partial c}{\partial z^{\prime }}+v^{\prime }\frac{\partial c}{\partial r^{\prime }}-\left(\frac{L}{U}\frac{\kappa }{L^2}\right)\frac{{\partial }^{2}c}{\partial z^{\prime 2}}.$$

$$\frac{1}{\epsilon }\mathcal{L}c=\frac{\partial c}{\partial t}+u\frac{\partial c}{\partial z}+v\frac{\partial c}{\partial r}-\epsilon \beta \frac{{\partial }^{2}c}{\partial z^2}.$$

$$\frac{\partial C}{\partial t}=G^{1}\Bigg{[}C,\frac{\partial C}{\partial z}\Bigg{]}+\varepsilon G^{2}\Bigg{[}C,\frac{\partial C}{\partial z},\frac{\partial^{2}C}{\partial z^{2}}\Bigg{]}+...=\sum_{n=1}^{\infty}\varepsilon^{n-1}G^{n}\Bigg{[}\left\{\frac{\partial^{j}C}{\partial z^{j}}\right\},j\in[0,n]\Bigg{]},$$

$$c=V^{0}\Bigg{[}r,t;C\Bigg{]}+\varepsilon V^{1}\Bigg{[}r,t;C,\frac{\partial C}{\partial z}\Bigg{]}+\varepsilon^{2}V^{2}\Bigg{[}r,t;C,\frac{\partial C}{\partial z},\frac{\partial^{2}C}{\partial z^{2}}\Bigg{]}+...=\sum_{n=0}^{\infty}\varepsilon^{n}V^{n}\left[r,t;\left\{\frac{\partial^{j}C}{\partial z^{j}}\right\},j\in[0,n]\right].$$

$$\begin{array}{rl}\frac{1}{\epsilon }\mathcal{L}{V}^0+\mathcal{L}{V}^1+\epsilon \mathcal{L}{V}^2\cdots =& \frac{d{V}^0}{dt}+u\frac{\partial V^0}{\partial z}+v\frac{\partial V^0}{\partial r}-\epsilon \beta \frac{{\partial }^2{V}^0}{\partial z^2}+\epsilon \frac{d{V}^1}{dt}+\epsilon u\frac{\partial V^1}{\partial z}+\epsilon v\frac{\partial V^1}{\partial r}-{\epsilon }^2\beta \frac{{\partial }^2{V}^1}{\partial z^2}+\dots \\ =& \frac{\partial V^0}{\partial t}+\frac{\partial V^0}{\partial C}(G^1+\epsilon G^2+...)+u\frac{\partial V^0}{\partial z}+v\frac{\partial V^0}{\partial r}-\epsilon \beta \frac{{\partial }^2{V}^0}{\partial z^2},\\ & +\epsilon \frac{\partial V^1}{\partial t}+\epsilon \frac{\partial V^1}{\partial C}(G^1+\epsilon G^2+...)+\epsilon \frac{\partial V^1}{\partial C^{(1)}}\frac{\partial (G^1+\epsilon G^2+...)}{\partial z},\\ & +\epsilon u\frac{\partial V^1}{\partial z}+\epsilon v\frac{\partial V^1}{\partial \eta }-{\epsilon }^2\beta \frac{{\partial }^2{V}^1}{\partial z^2}\dots \end{array}$$

$$\begin{array}{rl}\mathcal{L}{V}^0& =0,\dots \\ \mathcal{L}{V}^1& =\frac{\partial V^0}{\partial t}+\frac{\partial V^0}{\partial C}{G}^1+u\frac{\partial V^0}{\partial z}+v\frac{\partial V^0}{\partial r},\dots \\ \mathcal{L}{V}^n& =\frac{\partial V^{n-1}}{\partial t}+\sum _{\ell =1}^n\sum _{p=0}^{n-\ell }\frac{\partial V^{n-\ell }}{\partial C^{(p)}}\frac{{\partial }^p{G}^{\ell }}{\partial z^p}+u\frac{\partial V^{n-1}}{\partial z}+v\frac{\partial V^{n-1}}{\partial \eta }-\beta \frac{{\partial }^2{V}^{n-2}}{\partial z^2},\end{array}$$

$$\frac{\partial V^{n}}{\partial r}\bigg{|}_{r=1}=0\,\,\mathrm{for}\,\,\,\mathrm{all}\,\,n\geq 0.$$

$${\tau }_{\mathrm{adv}}=L/U\gtrsim a^2/\kappa ={\tau }_{\mathrm{diff}}^a$$

$$u=2U\left(1-\frac{r^2}{a^2}\right),$$

$$\kappa\frac{\partial c}{\partial r}\bigg{|}_{r=a}+\gamma c\big{|}_{r=a}=0.$$

$$\frac{\partial C}{\partial t}=-2\frac{\kappa }{a^2}(\gamma a/\kappa )\left(1-\frac{\gamma a/\kappa }{4}\right)C-\left(1+\frac{\gamma a/\kappa }{6}\right)U\frac{\partial C}{\partial z}+\kappa \left(1+\frac{a^2{U}^2}{48{\kappa }^2}\right)\frac{{\partial }^{2}C}{\partial z^2}.$$

$$\frac{\partial C}{\partial t}=$$

$$+\kappa \left(1+\frac{a^2{U}^2}{48{\kappa }^2}\left[1+\frac{(\gamma a/\kappa )}{30}\right]\right)\frac{{\partial }^{2}C}{\partial z^2}-\frac{a^4{U}^3}{2880{\kappa }^2}\frac{{\partial }^{3}C}{\partial z^3}.$$

$$\frac{\partial u}{\partial r}\bigg{|}_{r=a}=-\frac{a}{b}u(r=a),$$

$$u=2\alpha U_0\left(1-\frac{1}{\alpha }\frac{r^2}{a^2}\right),$$

$$\frac{\partial C}{\partial t}=-U\frac{\partial C}{\partial z}+\kappa \left(1+\frac{a^2{U}_0^2}{48{\kappa }^2}\right)\frac{{\partial }^{2}C}{\partial z^2}.$$

$$\frac{\partial C}{\partial t}=-U\frac{\partial C}{\partial z}+\kappa \left(1+\frac{a^2{U}^2}{48{\kappa }^2}I[u]\right)\frac{{\partial }^{2}C}{\partial z^2},$$

$$I[u]={\int }_0^{1}96r^{\prime }d{r}^{\prime }\left[\frac{u(r^{\prime })}{U}-1\right]\left[\frac{\tilde{u}(r^{\prime })}{U}-\frac{r^{\prime 2}}{4}\right],\tilde{u}(r^{\prime })={\int }_0^{r^{\prime }}dy{\int }_0^{y}dx\frac{xu(x)}{y}.$$

$$I[u]=\tilde{I}[\psi ]=96{\int }_0^1\frac{1}{x}{\left[{\int }_0^x{x}^{\prime }\psi (x^{\prime })d{x}^{\prime }\right]}^2.$$

$$u_{\mathrm{opt}}(r)=U\left(1+\alpha J_0\left(\frac{r}{a\sqrt{\lambda }}\right)\right).$$

$$u(r=a,t)=0,v(r=a,t)=\frac{\partial a}{\partial t}.$$

$$u(r,z,t)=2U\left(1-\frac{r^2}{a^2}\right),$$

$$v(r,z,t)=\frac{\partial a}{\partial t}\frac{r}{a}(2-\frac{r^2}{a^2})+2U\frac{\partial a}{\partial z}\frac{r}{a}(1-\frac{r^2}{a^2}),$$

$$U(z,t)=U(z=0,t)-\frac{1}{a^2}{\int }_0^{z}d\tilde{z}\frac{\partial a^2(\tilde{z},t)}{\partial t}.$$

$$\frac{\partial a^{2}C}{\partial t}=\frac{\partial }{\partial z}\left(-U{a}^{2}C+a^2\kappa \left(1+\frac{a^2{U}^2}{48{\kappa }^2}\right)\frac{\partial C}{\partial z}\right).$$

$$\frac{\partial \rho }{\partial t}=\frac{\partial }{\partial z}\left(-\left[U+\frac{\kappa }{a^2}\frac{\partial a^2}{\partial z}\left(1+\frac{a^2{U}^2}{48{\kappa }^2}\right)\right]\rho +\kappa \left(1+\frac{a^2{U}^2}{48{\kappa }^2}\right)\frac{\partial \rho }{\partial z}\right).$$

$$a(z,t{)}^2=a_0^2(1+\varphi \cos (\omega t-kz)).$$

$$a^2(z,t)U(z,t)=a^2(z=0,t)U(z=0,t)-\omega a_0^2\frac{\varphi }{k}[\cos (\omega t-kz)-\cos (\omega t)].$$

$$U=\frac{\varphi \omega }{k}\cos (\omega t-kz)\left(1-\varphi \cos (\omega t-kz)\right)+\mathcal{O}({\varphi }^2),$$

Full text

I Introduction

Microscale to nanoscale channels are fundamental building blocks for transport and separation of solutes in their stream. Interconnected channels form the vascular networks of animals, plants, fungi and slime molds - an essential transport system for all higher forms of life West (1999). Moreover, channels form lungs, intestines, the nephron of the kidney and the lymph system Freund et al. (2012); Cartwright et al. (2009). Channels are the key building block of microscale hydraulic engineering Stone et al. (2004), with particular importance for porous media applications Alim et al. (2017a) and in the emerging fields of smart materials Beebe et al. (2000) and soft robotics Shepherd et al. (2011). The dynamics of solute transport by flow in these channels is essential to transport signals Shapiro and Brenner (1988); Alim et al. (2017b), to power chemical reactions Zheng et al. (2017) and potentially program complex behavior Jubin et al. (2018). To increase the space of accessible behaviors one must go beyond the straight static version of a channel and incorporate the boundaries into design: carving for instance cross-sectionally varying channels Brenner (1980); Siwy (2006); Ajdari et al. (2006); Chinappi and Malgaretti (2018), or driving dynamics with active channel walls Marbach et al. (2018).

Yet, from the fluid physics perspective dispersion of solute in straight steady channels is already intricate. As solute particles inevitably transverse between streamlines of a standard Poiseuille flow profile by molecular diffusion, the effective dispersion of solute along the channel is enhanced. Taylor and Aris provided the basis to quantify this effect Taylor (1953, 1954); Aris (1956). The cross-sectionally averaged solute concentration $C(z,t)$ along the axis $z$ of a straight channel of radius $a$ with Poiseuille profile and cross-sectional averaged flow velocity $U$ can be approximated by

[TABLE]

where the bare molecular diffusion coefficient $\kappa$ is enhanced by a factor known as Taylor dispersion today, see Fig. 1. Since the Taylor dispersion effect is inevitable in microscale devices, significant effort has been invested to understand how to control and how to model Taylor dispersion within a number of settings.

Some of the key findings to extend Taylor dispersion in a number of situations include the understanding of the effect of reactive walls (a model for chemical reactions occurring at interfaces) Sankarasubramanian and Gill (1973); Lungu and Moffatt (1982); Smith (1983); Barton (1984); Shapiro and Brenner (1986); Balakotaiah et al. (1995); Biswas and Sen (2007); Datta and Ghosal (2008); Berezhkovskii and Skvortsov (2013); Levesque et al. (2013a, 2012, b), of buoyancy effects and cross flows (to model e.g. sedimenting particles) Erdogan and Chatwin (1967); Nadim et al. (1985); Shapiro and Brenner (1987); Lin and Shaqfeh (2019), of non-newtonian fluids Sankar et al. (2016) and extensions to more complex geometries such as porous media or networks Shapiro and Brenner (1988); Brenner (1980); Dorfman and Brenner (2002); Dutta and Leighton (2001); Dutta et al. (2006); Roberts (2004). In driven or more complex nonequilibrium systems a number of effects still remains to be explored (albeit the case of oscillatory shear flows as a model for the human blood vasculature Schmidt et al. (2005); Joshi et al. (1983)). For example, Taylor dispersion has entangled consequences in driven situations such as electro-osmosis or electro-phoretic transport Yariv and Dorfman (2007); Ng and Zhou (2012) and on the transport of active matter within porous media Alonso-Matilla et al. (2019); Bhattacharjee and Datta (2019). Yet, how one may control Taylor dispersion with active boundaries is still unanswered despite its utter importance for a number of biological systems, e.g. in situations that require mixing at very small scales Alim et al. (2017b); Cremer et al. (2016), or for fabrication of advanced soft robotics. We expect dispersion to be modified by complex flows at the inlet Ng and Zhou (2012) or dynamically pulsating boundaries Marbach et al. (2018), but to answer how, a thorough, quantitative understanding is required.

In this work we investigate how the effect of Taylor dispersion in straight channels can be diminished, suppressed or enhanced by specific boundary conditions: (1) controlling the flow profile entering the channel or (2) controlling the motion of the confining walls. To this end we employ the invariant manifold method to calculate the effective dispersion equation for the cross-sectionally averaged concentration profile Mercer and Roberts (1990, 1994). We show that changing the flow profile at the inlet can only serve to enhance dispersion, and the optimal flow profile for mixing is only a moderate increase as compared to classical Taylor dispersion in Poiseuille flow. However, we find that space-time contracting, active channel walls open up a broad range of possibilities to control Taylor dispersion, as they give rise to corrections additive to the Taylor dispersion correction. Thereby, we show that active, slow moving, channel walls allow to suppress Taylor dispersion and overall diminish dispersion. Very fast walls allow to enhance Taylor dispersion by shuttle dispersion. We fully rationalize all these regimes and their crossover.

Our paper is organized in the following way. In the perspective of making a pedagogical introduction to the invariant manifold method used to derive asymptotic dispersion equations for the solute concentration profile, we dedicate section II to the explanation of the method and to an example. We revisit this rigorous method and show how versatile and powerful it can be, as it is easily applied to compute analytic results for the various settings that we explore here. Notably, we stress that while the standard textbook approach to derive Taylor dispersion may fail in more advanced settings, the invariant manifold expansion robustly yields reliable results. Section III reports the results on the impact of flow profile in straight channels, notably asking the question of how slippage modifies dispersion. Section IV reports the results on the impact of channel wall space-time contractions. We report in the Appendix the detailed calculations for each section.

II Background: the invariant manifold method

Dispersion of solute of concentration $c(r,z,t)$ in a long, slender channel of radius $a(z,t)$ and length $L$ is fully described by the advection-diffusion equation

[TABLE]

where $u(r,z,t)$ and $v(r,z,t)$ are the axial and radial flow velocities and $\kappa$ is the molecular diffusivity. To solve Eq. (2) one must also specify boundary conditions for the solute flux at the channel wall. While the advection-diffusion equation is the most accurate description of dispersion, the impact that flow profiles or boundary conditions have on the effective transport and dispersion is hard to infer as two spatial components, $r$ and $z$, have to be considered. The very insightful approach by Taylor and Aris Taylor (1953, 1954); Aris (1956) has been to instead aim for a reformulation of the problem in terms of the dispersion of the cross-sectionally averaged solute concentration $C(z,t)=\frac{2}{R^{2}}\int c(r,z,t)\,r\,dr$, averaging out the dynamics in radial direction. While heuristically averaging out the radial direction works for dispersion in a straight channel with Poiseuille flow profile (see Appendix A) more complicated flow profiles, e.g. dependent in $z$ and $t$, require more formal techniques.

Among those techniques the invariant manifold method originally introduced by Mercer and Roberts Mercer and Roberts (1994, 1990) is very versatile and can be readily applied to derive dispersion dynamics in long slender channels for any boundary condition and flow profile, as outlined in the following. Notably, the invariant manifold method was recently used as a rigorous proof to establish the effective diffusion expression in a straight channel Beck et al. (2018). The method is based on an invariant manifold description, which has the advantage to be systematically extendable to higher orders, and thus allows to keep faithfully track of the appropriate order required and also allows to assess the magnitude of the next higher order term neglected. For more information on the invariant manifold description and other methods we refer the reader to Refs. Carr (1981); Roberts (1989); Balakotaiah et al. (1995); Watt and Roberts (1995); Rosencrans (1997).

Although there exists a number of other techniques to infer the effective diffusion coefficient at long times (for example the method of moments, implicit constructions and others Gill and Sankarasubramanian (1970); Roberts (2014); Frankel and Brenner (1989)), the invariant manifold method is one of the only ones that can provide asymptotic dispersion equations (and not just dispersion coefficients), which is one of our aims. We refer the reader to Ref. Young and Jones (1991) for a detailed comparison of the different methods.

II.1 Overview of the invariant manifold method

As a first step to introduce the invariant manifold method we rewrite Eq. (2) by introducing non-dimensional quantities: $r^{\prime}=r/a$ as the radial coordinate, $z^{\prime}=z/L$ the axial coordinate, $t^{\prime}=t/(L/U)$ time with $U$ the mean longitudinal flow velocity in the channel, $u^{\prime}=u/U$ as the longitudinal velocity, and $v^{\prime}=v/(Ua/L)$ the radial velocity. The advection-diffusion equation to be expanded is then given, after reordering terms, by:

[TABLE]

In further equations we will drop the prime notations for simplicity. In the following we also use the abbreviated operator $\mathcal{L}c=\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial c}{\partial r}\right)$.

The aim is now to find an asymptotic dispersion equation. This is done typically by assuming that dynamics of solute dispersion across the channel are much faster than along the channel. More formally this means that the time to diffuse across the cross section $\tau_{\rm diff}^{a}=a^{2}/\kappa$ is much smaller than the time to be advected along the channel $\tau_{\rm adv}=L/U$. We can therefore introduce a small parameter $\varepsilon=\tau_{\rm diff}^{a}/\tau_{\rm adv}$ and look for an expansion in this small parameter 111It has further been shown that the invariant manifold method can hold even in the case where this small parameter approaches unity Rosencrans (1997).. Furthermore, we assume that the time to diffuse along the channel main axis $\tau_{\rm diff}^{L}=L^{2}/\kappa$ is much longer than $\tau_{\rm adv}$ and that we may write $\tau_{\rm diff}^{L}=\beta\varepsilon$ where $\beta$ is of order unity. The advection-diffusion equation simplifies to:

[TABLE]

We now look for an effective equation on the cross-sectionally averaged concentration $C(z,t)$ and we aim for a systematic expansion in $\varepsilon$,

[TABLE]

where the $G^{n}$ represent the different terms in the asymptotic differential equation when it is extended to order $n$ in $\varepsilon$. Note that because of the structure of the differential equation, $G^{n}$ will only depend on the spatial derivatives of $C$ in $z$ of order $n$ or less. For the solute concentration we also look for a systematic expansion as

[TABLE]

To infer the solution for each successive order of $V^{n}$ and $G^{n}$ we first substitute both Eqs. (5) and (6) into Eq. (4). Note, that we have to be particularly careful in respecting the chain rule when taking the time derivative of $V^{n}$ since $V^{n}$ not only directly depends on time but also indirectly via it’s dependence on all individual $\frac{\partial^{j}C(z,t)}{\partial z^{j}}$:

[TABLE]

Now we may equate terms of the same order in $\varepsilon$ and obtain:

[TABLE]

where $C^{(p)}=\frac{\partial^{p}C}{\partial z^{p}}$. Note that $v\frac{\partial V^{n-1}}{\partial\eta}$ is of the same order as $u\frac{\partial V^{n-1}}{\partial z}$, this can be seen from the fact that the divergence of the flow field is zero.

In addition $V^{n}$ has to obey a boundary condition at the channel wall. In the simplest case, we expect an impermeable wall and therefore the concentration profile has to obey a no flux boundary condition \frac{\partial c}{\partial r}\big{|}_{r=1}=0. This translates into

[TABLE]

Integration constants when solving the equations are further constraint as $V^{n}$ and its radial derivatives need to stay finite at $r=0$. Last, a simple closing equation has to be verified, namely $\int 2V^{0}r\,dr=C$ and $\int V^{n}r\,dr=0$ for all $n\geq 1$.

Solving at each order Eq. (8) and replacing the $G^{n}$ terms in Eq. (5) yields the desired asymptotic dispersion equation. Finally, it is of interest to truncate the expansion at finite order (most commonly at $n=2$) to be able to use a simple analytic description of dispersion. This truncation is justified since the asymptotic series is convergent Mercer and Roberts (1990, 1994); Beck et al. (2018) provided that $\varepsilon\gtrsim 1$, which corresponds to the following physical condition

[TABLE]

where $L$ is the typical length scale (along $z$) of observation, $a$ the channel radius, $\kappa$ the molecular diffusivity and $U$ the typical longitudinal velocity Mercer and Roberts (1994). In the following we will revert to fully dimensional equations.

Employing the invariant manifold method on the classical example of Taylor dispersion in a straight channel with steady Poiseuille flow profile yields at second order (in $\varepsilon$) the classical result of Eq. 1 (see Appendix B). The same result can be derived in a more heuristic way (see Appendix A). As stated earlier, more complicated dispersion dynamics are not well captured by the heuristic approach as exemplified in the following subsection.

II.2 Example: impact of absorption at the channel wall on dispersion in a straight channel

The power of the invariant manifold approach versus the heuristic approach is best exemplified for a straight channel with solute absorption at the channel wall. The problem has recently been solved with the heuristic approach Meigel and Alim (2018), at the same time a number of theoretical expansions and numerical solutions are available Sankarasubramanian and Gill (1973); Lungu and Moffatt (1982); Smith (1983); Barton (1984); Shapiro and Brenner (1986); Balakotaiah et al. (1995); Biswas and Sen (2007); Datta and Ghosal (2008); Berezhkovskii and Skvortsov (2013); Levesque et al. (2013a, 2012, b). Among them we use Ref. Balakotaiah et al. (1995) as a benchmark for the exact expansion coefficients.

The fluid is flowing in a Poiseuille profile down the channel, with zero radial flow component, $v=0$, and longitudinal flow velocity

[TABLE]

where $U$ denotes the cross-sectionally averaged flow velocity. Solute is absorbed passively at the channel wall, see Fig. 2:

[TABLE]

where $\gamma$ is the surface absorption coefficient in units of $\mathrm{m.\,s^{-1}}$. To solve step by step the expansion coefficients with the invariant manifold approach it is important to define how the different orders of $V^{n}$ are affected by the solute flux boundary condition. Since we are interested in the cross-sectionally averaged concentration $C(z,t)$ as a proxi for $c(r,z,t)$, we restrict ourselves to the case where cross-sectional gradients in absorption are averaged out quickly by cross-sectional diffusion of solute,* i.e.* $\frac{a}{\gamma}\gg\frac{a^{2}}{\kappa}$, such that $\gamma a/\kappa\ll 1$. Furthermore, we assume that $\gamma a/\kappa$ is of the same order as the small number $\varepsilon$ describing the expansion of the invariant manifold method. As a consequence, we define the boundary conditions as \kappa\frac{\partial V^{n}}{\partial r}\big{|}_{r=a}=-\gamma V^{n-1}|_{r=a} for $n\geq 1$ and \frac{\partial V^{0}}{\partial r}\big{|}_{r=a}=0. Working out the expansion to second order, see Appendix C, we find,

[TABLE]

In Eq. (11), the first term on the right hand side is a sink term corresponding to the absorption at the channel wall. The expansion allows to show that absorption increases the effective drift at second order by $U\left(1+\frac{\gamma a/\kappa}{6}\right)$ . In fact, as particles being absorbed are closer to the walls, the particles’ center of mass is effectively pushed downstream. The effective diffusion constant is exactly that of Taylor’s dispersion in a straight channel, so absorption at second order does not impact diffusion here.

Beyond the mechanistic understanding, it is insightful to compare Eq. (11) to the equation resulting from the heuristic approach Meigel and Alim (2018) (or see Appendix C) and the numerical exact solutions Balakotaiah et al. (1995), see Fig. 3.

We find that the heuristic approach fails at predicting the drift and the effective diffusion accurately for moderately small absorption rates, see Fig. 3. Overall the second order invariant manifold expansion, Eq. (11), only starts to show deviations from the exact results at large absorption parameters. This allows us to conclude that any heuristic derivation of slightly complex Taylor dispersion problems can be made in an unreliable way, and that more advanced techniques such as the method of moments or the invariant manifold method should be used instead.

Note that one may push the expansion to third order, see Appendix C, to capture corrections in particular to the effective diffusion,

[TABLE]

Here, we see that if absorption is small, dispersion can be enhanced by absorption. Since absorption reduces solute concentration at the channel wall, it helps to increase the radial gradient of solute concentration and thus contributes to enhancing the redistribution effect. At even higher absorption rates – pushing the expansion to fourth order – solute concentration is reduced near the wall too rapidly for any redistribution of solute between slow and fast streamlines to take place, thus diminishing Taylor dispersion again. The key insight gained here is that Taylor dispersion can be enhanced by changing the concentration profile of the solute across the channel’s cross-section with e.g. absorbing boundary conditions - although the magnitude of the effect is relatively small compared to pure Taylor dispersion. To investigate other means to control Taylor dispersion we now turn to the role of the flow profile in Taylor dispersion.

III Result 1: Impact of flow profile in a straight channel

on Taylor dispersion

The question of how to prescribe the inlet flow to control the downstream dispersion of solutes is of general interest for a number of industrial applications, and is part of the broader question of controllability of flows Glass and Horsin (2012, 2016). We here investigate the impact of the axial flow profile on Taylor dispersion. We will therefore deviate from pure Poiseuille profile by first investigating the impact of slip boundaries and subsequently deriving the optimal axial flow profile for enhancing Taylor dispersion.

III.1 Impact of slip-boundary condition at the channel wall

In the context of microfluidics, it is interesting to focus on flow profiles that take into account different boundary conditions. In fact in micro to nano systems, the boundary properties may largely influence the flow profile in the bulk because the surface to volume ratio increases. For instance, electro-osmotic or diffusio-osmotic flows may occur, and solute concentration may also produce a feedback on the flow Tabeling (2009). A standard scenario arising in microfluidics Lauga et al. (2007) and with even more impact in nanofluidics Secchi et al. (2016) is slip at the boundary of the channel wall.

Imposing the Navier slip boundary condition,

[TABLE]

where $b>0$ is a length called the slip length, results in a vanishing radial flow velocity $v\equiv 0$ and an axial flow velocity in the channel

[TABLE]

where $\alpha=1+2b/a$. $U_{0}$ is the cross-sectionally averaged axial flow velocity field in the absence of slip, characterizing a constant pressure drop through the Poiseuille law U_{0}=\frac{a^{2}}{8\eta}\big{|}\frac{\partial p}{\partial z}\big{|}. Note, that the acting cross-sectionally averaged axial flow velocity $U=U_{0}(2\alpha-1)=U_{0}\left(1+\frac{4b}{a}\right)$ is bigger than the flow velocity without slip.

Solving for Taylor dispersion with the invariant manifold expansion to second order, see Appendix D-1, we find,

[TABLE]

Eq. (14) converges to the expected solution with no slip at the boundary Eq. (40), when $b\rightarrow 0$, as required. The effective drift that pushes the solute is indeed $U$ and not $U_{0}$, yet the effective diffusion coefficient is independent of the slip length. The part that counts for enhancement of dispersion is exactly the deviation of the flow profile from a plug flow. In fact, the flow velocity deviation to plug flow, $u(r)-\min_{r}u(r)=2U_{0}(1-r^{2}/a^{2})$, corresponds to the associated Poiseuille flow with no slip at the channel wall and thus gives exactly the same dispersion enhancement as standard Poiseuille flow. In more complex situations however, coupled effects may induce a positive feedback of slippage on Taylor dispersion, e.g. for Taylor dispersion by electro-osmotic flow that is itself increased by slippage Ng and Zhou (2012).

III.2 Impact of an arbitrary axial flow profile

To investigate if any axial flow profile can outcompete Taylor dispersion as arising from a Poiseuille flow profile we turn to derive Taylor dispersion for an arbitrary flow profile. This is useful to model transport through complex fluids such as Herschel–Bulkley fluids, or to perform advanced design of devices. For example one could tailor the inlet flow profile in a microfluidic device with an array of inlets with different fluid flows Lin et al. (2004). Here, the advantage of the invariant manifold method is important, in that the method is applicable independent of the characteristics of the flow profile, and thus may be used for any flow profile. Compared to similar approaches Beck et al. (2018) we find an analytic expression for the Taylor dispersion coefficient, readily applicable to any flow profile.

For a given axial flow profile $u(r)$ and zero radial flow velocity $v=0$ we find, see Appendix D-2:

[TABLE]

where $U=2\int u(r^{\prime})r^{\prime}\,dr^{\prime}$ denotes the cross-sectionally averaged flow velocity, with $r^{\prime}=r/a$ the non-dimensional radial coordinate, and $I[u]$ is a correction factor to Taylor dispersion that is defined by

[TABLE]

As expected, the solute is effectively advected by $U$, the cross-sectionally averaged flow velocity over the cross-section. Moreover enhancement to the effective dispersion is present via the Taylor mechanism. $I[u]$ is a measure of the strength of this correction. In this formulation it is straight-forward to see that Taylor dispersion enhancement vanishes for plug flow, e.g. if $u(r^{\prime})\equiv U$ then immediately one finds $I[u]=0$. If the profile is Poiseuille-like then $I[u]=1$. It is insightful to rewrite the correction factor by introducing the non-dimensional function $\psi(r^{\prime})=\frac{u(r^{\prime})}{U}-1$. As $\int_{0}^{1}r^{\prime}\,dr^{\prime}\,\psi(r^{\prime})=0$ this function measures the deviation of the velocity profile $u$ to a plug flow. It allows to rewrite the correction factor $I[u]$ as

[TABLE]

One immediately sees that the integrand is always positive and thus there can never be solute focusing by flow in a straight channel. Any flow profile that deviates from plug flow results in $I[u]>0$ an thus gives rise to enhanced dispersion. Yet, the remaining question is how big $I[u]$ can get.

III.3 Optimal axial flow profile to enhance Taylor dispersion

We now use the previous results (Eqs. (15) and (16)) to search for the optimal flow profile that maximizes $\tilde{I}[\psi]$ without changing the average flow velocity. Hence, we assume that the cross-sectional average of $\psi(r^{\prime})$ is zero and that the deviation from plug flow has the same strength as Poiseuille flow, i.e. the same $L^{2}$ norm. Under these constraints we derive the optimal axial flow, see Appendix D-3, to be given by a Bessel function $\psi(r^{\prime})=\alpha J_{0}\left(\frac{r^{\prime}}{\sqrt{\lambda}}\right)$ with $\lambda\simeq 0.068$ and $\alpha\simeq 1.43$. Translated back to axial flow velocities this implies

[TABLE]

The optimal profile is much steeper than the Poiseuille flow, particularly with two different regions, see Fig. 5: the central region moves really fast and the borders move slowly. This gradient will maximize the transverse diffusive drift and increase the effective lateral dispersion. That said, $I[u_{\rm opt}]\simeq 1.09$ is only slightly larger than $I[u]=1$ obtained for Poiseuille flow. Obviously, one would have to work with a very specific setting or a very peculiar fluid to obtain such a flow profile. Yet, it is fascinating to see that Poiseuille flow is almost the optimal flow profile for Taylor dispersion. In the same context it seems insightful to evaluate also the correction factor for the most common description of blood, namely a Herschel-Bulkley fluid. With $r_{p}=0.2$ and $n=0.95$ (taking $r_{p}$ and $n$ as parameters of the Herschel-Bulkley fluid from Sankar et al. (2016)) the Taylor dispersion enhancement factor is $I[u_{\mathrm{H-B}}]\simeq 0.96$ and with $n=1.05$, $I[u_{\mathrm{H-B}}]\simeq 0.91$. Here, the flat center of the Herschel-Bulkley flow profile prevents the Taylor-dispersion effect in the center, lowering the overall enhancement of dispersion Sankar et al. (2016). While the flow profile clearly controls Taylor dispersion, enhancing Taylor dispersion merely on the basis of the axial flow profile is limited. We, therefore, next turn to investigate the impact of active channel walls.

IV Result 2: Impact of channel wall space-time contractions on Taylor dispersion

Turning to active channel walls we here aim to derive Taylor dispersion in general for any space-time contraction dynamics of a channel radius $a(z,t)$ going beyond results for spatially varying channels Mercer and Roberts (1994) and the role of oscillatory shuttle flow Chatwin (1975); Watson (1983); Leighton and McCready (1988); Schmidt et al. (2005); Joshi et al. (1983). The channel contractions now imply a radial component of fluid flow due to the no-slip boundary condition at the channel wall imposing

[TABLE]

Taking into account the continuity equation $\frac{\partial u}{\partial z}=-\frac{1}{r}\frac{\partial}{\partial r}\left(rv\right)$ fully defines the flow velocities

[TABLE]

where the cross-sectionally averaged axial flow velocity $U$ follows from conservation of mass $\frac{\partial}{\partial z}\left(a^{2}U\right)=-\frac{\partial}{\partial t}a^{2}$, resulting in

[TABLE]

IV.1 Long time concentration evolution in a space-time contracting channel

Using the invariant manifold method up to second order we find for the dynamics of the cross-sectionally averaged concentration $C$ (see Appendix E),

[TABLE]

Eq. (20) is eventually remarkably simple and resembles the classical Taylor Dispersion equation in the case of a steady straight channel Eq. (1). With the notation of Eq. (1) we may identify as the relevant variable, $\rho=a^{2}C$ the particle mass in an infinitesimal slice (equivalent to the marginal probability to find solute particles at coordinate $z$ and $t$). In Eq. (20) particle mass is advected with the effective drift $U$, the cross-sectionally averaged axial flow velocity, while effective diffusion is proportional to the gradient of the solute itself $\frac{\partial C}{\partial z}$. One may readily use Eq. (20) to model solute dispersion in dynamic channels (or networks of channels), not limited to but particularly prevalent, in living systems.

Rewriting one more time in terms of particle mass $\rho$ only, we find

[TABLE]

In this formulation we see that the particle mass’s drift is modulated by the deformation of the channel $\frac{\partial a^{2}}{\partial z}$. For a pure Brownian particle in a channel in the absence of any flow the impact of confinement shape on a particle’s effective’s drift is known as entropic slow down Martens et al. (2013). Here, we find that flow velocity modulates entropic slow down and may augment the effect. Yet, how strongly the Taylor dispersion term $\frac{a^{2}U^{2}}{48\kappa}$ changes effective drift and effective dispersion cannot be estimated in this very general form, since here both $a(z,t)$ and $U(z,t)$ depend on space and time. We therefore turn to the example of peristaltic pumping to estimate the impact of all terms on effective dispersion.

IV.2 The impact of peristaltic pumping on dispersion with emphasis on the role of Taylor dispersion

To estimate the impact of peristaltic pumping on Taylor dispersion we now consider explicit space-time contraction of the channel wall of the functional form

[TABLE]

Note, that to first order in $\phi\ll 1$ this ansatz is equivalent to the classical example of peristaltic pumping Shapiro et al. (1969), yet this form is analytically more tractable Alim et al. (2013). From conservation of mass, Eq. (19), it follows that

[TABLE]

To allow a fully analytic calculation to gain insight into the process, we need to simplify the space-time dependence of the mean flow $U$. To this end, we consider an oscillatory inflow boundary condition $Q(z=0,t)=a^{2}U(z=0,t)=-\omega a_{0}^{2}\frac{\phi}{k}\cos(\omega t)$ accounting for the flow generated elsewhere in a long tube. With this assumption we find

[TABLE]

where we used the fact that $\phi\ll 1$. We further shift our coordinate system by performing the change of variables: $\xi\leftarrow\frac{\omega}{k}t-z$ and $t\leftarrow t$ such that Eq. (21) becomes a dissociable function of $\xi$ and $t$ as

[TABLE]

where now advection and drift terms all depend only on $\xi$. Eq. (25) is now amenable to further analytic computations. Either based on the calculation of mean first passage times Reimann et al. (2001) or on the analysis of probability distribution functions Guérin and Dean (2015) it is possible to derive the effective dispersion coefficient at infinitely long times, see Appendix F, arriving at,

[TABLE]

where we here adopted the Taylor dispersion notation in terms of the Péclet number. To allow comparison to the case of the straight steady channel, we here define a time averaged Péclet number as $Pe=\sqrt{\frac{\langle U^{2}a^{2}\rangle_{T}}{\kappa^{2}}}=\frac{\phi\omega a_{0}}{\sqrt{2}k\kappa}$, where $\langle.\rangle_{T}$ indicates the time-average over a contraction period $T=2\pi/\omega$. We identify two correction terms to the Taylor dispersion term in $Pe^{2}/48$.

The first correction, $\frac{1}{2}\phi^{2}\frac{6Pe^{2}}{2Pe^{2}+(k\phi a_{0})^{2}}$, is a shear dispersion contribution. It arises as the alternating shear pushes particles back and forth, and as they diffuse on top of that, effectively increasing effective dispersion. This effect was identified separately in the case of periodic shuttle flows in straight channels Chatwin (1975); Watson (1983); Leighton and McCready (1988); Schmidt et al. (2005); Joshi et al. (1983) and is sometimes referred to as shear dispersion just like Taylor dispersion. We must stress that the two effects are entirely different: Taylor dispersion arises because of shear in the lateral direction, so because of a non-uniform flow profile in the lateral direction but arises independently of dispersion enhanced by shuttle streaming. We therefore refer here to the later effect as shuttle dispersion.

The second correction, $-\frac{1}{2}\phi^{2}\frac{(k\phi a_{0})^{2}}{2Pe^{2}+(k\phi a_{0})^{2}}$ is a negative contribution to effective dispersion. It occurs as the channel walls are oscillating in space, therefore creating constrictions that are hard to transverse and cavities that act as particle traps. This effect counteracts dispersion, a phenomenon called entropic slow down. It was previously identified in the absence of temporal variations, or more formally when $Pe=0$, by numerous authors Burada et al. (2009); Martens et al. (2013); Malgaretti et al. (2013); Yang et al. (2017) and is here generalized for additional flow. Remarkably, all these contributions to effective dispersion are additive, which is not necessarily obvious. This is particularly striking since all of these effects are different in their origin and may arise independently. They are now easily comparable through Eq. (26).

To estimate the magnitude of the correction terms relative to the Taylor dispersion coefficient in a straight channel it is insightful to rewrite the effective dispersion in terms of a Péclet number for advective versus diffusive transport along the contraction wave of wave number $k$, $Pe_{\mathrm{wave}}=Pe/(k\phi a_{0})=\omega/(\sqrt{2}\kappa k^{2})$,

[TABLE]

In this notation we see that all corrections to molecular diffusivity $\kappa$ scale with $\phi^{2}$, yet the Taylor dispersion correction has an additional prefactor $(ka_{0})^{2}$ which is always smaller than one and particularly small for long wavelengths of the space-time contractions. Exploring the limits, we find that for $Pe_{\mathrm{wave}}\ll 1$ – a slowly contracting channel –, shuttle dispersion and Taylor dispersion corrections vanish, with Taylor dispersion vanishing faster. Entropic slow down dominates and simplifies to $-\frac{1}{2}\phi^{2}$. For more intermediate values of $Pe_{\mathrm{wave}}$ entropic slow down vanishes and shuttle dispersion simplifies to $\frac{3}{2}\phi^{2}$, while Taylor dispersion may not yet be relevant. For very large $Pe_{\mathrm{wave}}$ however, the Taylor dispersion correction again dominates and continues to increase with $Pe_{\mathrm{wave}}$ as $Pe_{\mathrm{wave}}^{2}$. The different regimes are illustrated in Fig. 7, that clearly demonstrates that different parameters will yield very different corrections to dispersion. Although the transition between entropic slow down and shuttle dispersion was clearly evidenced by one of the authors Marbach et al. (2018), it is now striking to see how Taylor dispersion dominates either regime in a range of parameters, typically when $ka_{0}\gtrsim 8/Pe_{\rm wave}$ or more simply when $Pe\gtrsim 8$. Note, this requirement is also compatible with the validity condition of the framework Eq. (10) as long as $L/a_{0}\gg 1$. The results displayed in Fig. 7 suggest a number of ways to control and tune dispersion in an active channel.

With such peristaltic pumping, one finds also that the fluid within the channel is advected with a mean non-zero effective drift (see Appendix F) that writes

[TABLE]

When the contractions disappear, naturally this effective drift disappears. Making the contractions bigger (i.e. $\phi$ bigger) increases the effective drift (by increasing the displaced mass). The effective drift is in the direction of the peristaltic wave and scales with the phase velocity $\omega/k$. One may therefore control the direction and the amplitude of the effective drift by tuning the parameters of the peristaltic pump. One may also keep the effective drift constant while changing the effective dispersion (see Fig. 7-b). Note that one may change further this effective velocity by adapting the boundary conditions at the inlet and outlet of the tube.

V Discussion and Conclusion

The recent shift to active and smart artificial systems requires to revisit Taylor dispersion afresh and explore the impact of specific boundary conditions on the walls and extend to active channel walls. Moreover, we now have the ability to design innovative systems and to make refined observations with emerging techniques that allow us to quantify flows at very small scales in various biological contexts – flows generated by acto-myosin cortical activity in amoeba Lewis et al. (2015), slime molds Matsumoto et al. (2008) or plant cells Peremyslov et al. (2015) or flows generated by cilia motion Jonas et al. (2011) – but also in nanoscale artificial systems down to a few fL/s Secchi et al. (2016). Therefore, the ability to model how specific flows and boundary conditions affect solute transport within channels is critical to understand solute transport through these systems.

By revisiting a rigorous and systematic method, the invariant manifold approach, we showed that we could assess the impact of various boundary conditions (solute absorption, flow slippage), flow profile and active channel walls on long time advection and dispersion of the solute. We stress here again that the invariant manifold approach provides a rigorous framework that allows to derive an effective advection-diffusion equation for the solute, along the principal axis of interest (the long axis of a tube for example) in a systematic way. Importantly, this method avoids averaging mistakes in more complex situations than just a straight tube (as for absorption boundary conditions for example), and also allows to build more refined approximations with e.g. higher order derivatives in the effective advection-diffusion equation.

We recapitulate the specific results of our work:

(1) controlling the flow profile entering the channel gives a small handle to tune the effective dispersion. In particular, we have shown that allowing for slip on the channel walls, but still with a constant pressure drop along the channel does not modify dispersion. Furthermore, the Poiseuille flow profile in a cylindrical channel is nearly optimal for dispersion. As flow profiles get “flatter”, the dispersion enhancement diminishes to zero (zero being obtained for a plug flow profile).

(2) controlling the motion of the confining walls gives a huge handle to tune both effective dispersion and effective drift of the solute. The average solute concentration profile $C(z,t)$ at location $z$ on the principal axis of the tube may be described by

[TABLE]

where $a(z,t)$ is the radius of the pulsating confining wall and $U$ is the mean fluid velocity triggered by the contractions and flow incompressibility. This effective model for the concentration profile evolution is valid as long as Eq. (10) is verified, that we recall here for consistency: $\frac{L}{U}\gtrsim\frac{a^{2}}{\kappa}$. Note that it is possible to go beyond this regime, by pursuing the expansion and allowing for higher order derivatives in Eq. 29. Furthermore, for active channel walls, we showed that it is possible to tune dispersion between three different regimes of dispersion arising in a parameter space with two degrees of freedom: the contraction frequency and the contraction wavelength of the active channel walls.

Our work highlights the clear versatility of the invariant manifold method. It would be insightful to see the method be used to investigate further settings beyond the scope of this work for example coupling various effects, and to obtain also simple and readily applicable results for instance to pinpoint the role of diffusio-osmotic or electro-osmotic flows on dispersion Ng and Zhou (2012) with different boundary conditions.

Beyond the regimes that we have considered – where our continuum model is applicable – a number of extensions are possible. For example, it would be enlightening to understand how large peristaltic contractions affect the effective dispersion (with $\ell\simeq a$ where $\ell$ is the typical wavelength of the contraction). This could be accessible via continuum simulations with realistic membrane models for the tubes Deserno (2015). Furthermore, when the system size drops to nanoscales, continuum theories will fail to describe in a reliable way the velocity field in the fluid, and other effects. Only brownian dynamics or molecular dynamics are able to access these scales, assess the behaviors at play Yoshida et al. (2018) and inform us on the possible mechanisms by which dispersion may be enhanced or decreased. More efforts are needed in particular to bridge the gap between those discrete scales and continuum scales.

Acknowledgements.

The authors are indebted to Lydéric Bocquet, Agnese Codutti, Michael P. Brenner, David S. Dean, Frédéric Marbach and Felix J. Meigel for fruitful discussions. S. M. acknowleges funding from A.N.R. Neptune. K. A. acknowledges funding from the Max Planck Society.

Appendix

V.1 Textbook derivation of Taylor dispersion

Here we revisit for completeness Taylor dispersion as it is typically introduced in a textbook (e.g. Bruus (2008)). Consider solute dispersing in a straight cylindrical channel of radius $a$ with radial and longitudinal coordinates $(r,z)$, see Fig. 8.

At low Reynolds numbers, the flow field induced by a pressure drop in the channel obeys a Poiseuille profile

[TABLE]

where $U$ is the cross-sectional average of the flow velocity $u(r)$. Apart from this definition of $U$ we will use an overline $\overline{f}$ to denote the cross-sectional average of a variable $f$ in the following. Suppose at time $t=0$, an axisymmetric solute concentration is released $c(z,r,t)$. Then the solute concentration obeys the following advection diffusion equation:

[TABLE]

We look for an equation for $C(z,t)=\overline{c(z,r,t)}$ assuming that $\tilde{c}(z,r,t)=c(z,r,t)-C(z,t)$ is a small correction to the cross-sectional average concentration. For coherence with the following section we also write $C(z,t)=C(z,t)$. Inserting that in Eq. (30) gives

[TABLE]

and taking the cross-sectional average and using the boundary conditions gives

[TABLE]

Eq. (32) is the partial differential equation for $C$ that we are looking for. The problem is to evaluate $\overline{u(r)\frac{\partial\tilde{c}}{\partial z}}$. Subtracting Eq. (32) from Eq. (31) gives:

[TABLE]

We can use this expression to solve for the missing term but we have to perform coarse grain approximations. First we look for “long time” solutions, typically for times longer than the typical time to diffuse across the cross-section $t\gg\frac{a^{2}}{\kappa}$, where we expect that the cross channel diffusion will have averaged out such that we may neglect $\tilde{c}_{t}$. For such time scales we also expect $\tilde{c}\ll C$, such that we may also neglect $\left(u(r)\frac{\partial\tilde{c}}{\partial z}-\overline{u(r)\frac{\partial\tilde{c}}{\partial z}}\right)$. Furthermore, since $\tilde{c}$ is essentially a radial correction we assume that gradients in the $r$-direction for $\tilde{c}$ will be much greater than those in the $z$-direction. As a result Eq. (33) reduces to:

[TABLE]

Integrating, using that $\tilde{c}$ is finite at $r=0$ and the fact that the cross-sectional average of $\tilde{c}$ should vanish yields:

[TABLE]

Then one is left to perform the cross-section average of $u(r)\frac{\partial\tilde{c}}{\partial z}$ and one finds rather easily that Eq. (32) can be written as:

[TABLE]

This equation is called the Taylor dispersion equation. We notice that it is a simple advection diffusion equation along the axial direction, with an effective drift that is exactly the average drift, and an effective diffusion that is enhanced as compared to the bare diffusion $\kappa$. This enhancement factor is due to the fact that the Poiseuille flow is not a plug flow, and thus advects particles at different radial positions with different speeds. After a typical time $\tau$ the displacement will be $\delta\ell=U\tau$. This induces a radial and axial shift in the dispersion profile, that generates radial diffusion gradients. To equilibrate over the cross section, the time $\tau$ has to match the time $a^{2}/\kappa$. As a result the effective diffusion due to the Poiseuille flow is $\delta\ell^{2}/\tau=\frac{U^{2}a^{2}}{\kappa}$, which is exactly the enhancement factor deduced analytically in Eq. (36).

V.2 Invariant manifold method for classical Taylor dispersion in a straight channel

We here derive Taylor dispersion again for a straight channel with Poiseuille flow profile, but now with the invariant manifold method. $U$ stands for the cross-sectional average of $u$, and the fluid flow reads:

[TABLE]

where now $\eta=r/a$ is the non-dimensional radial coordinate (to lighten notations compared to $r^{\prime}$ in the main text), and the only non-dimensional variable that we keep for the derivation.

0th order expansion

We first solve for $V^{0}$. $\mathcal{L}V^{0}=0$ with \frac{\partial V^{0}}{\partial\eta}\big{|}_{\eta=1}=0 gives $V^{0}=C$. Note, that $V^{0}$ does not depend explicitly on $t$ nor $\eta$, therefore $\frac{\partial V^{0}}{\partial t}=0$.

1st order expansion

Now we solve for $V^{1}$. The defining equations are:

[TABLE]

We use the expression for $u(\eta)$ and integrate once with respect to $\eta$:

[TABLE]

where $A$ is an integration constant. It vanishes since $V^{1}$ and its derivative have to stay finite at $\eta=0$. Using the boundary condition at $\eta=1$, yields the expression for $G^{1}$:

[TABLE]

We now use this to derive an expression for $V^{1}$, and integrate once more. The integration constant is found using the fact that the average over the cross section of $V^{1}$ is zero, $\int V^{1}\eta\,d\eta=0$, such that :

[TABLE]

2nd order expansion

Now we solve for $V^{2}$. The defining equations are:

[TABLE]

Now we use Eq. (38) and the expression for $u$; and integrate once:

[TABLE]

where the integration constant vanishes since $V^{2}$ has to be finite at $\eta=0$, similarly to the previous integration. Using the boundary condition for $V^{1}$ yields:

[TABLE]

Combining this result with Eq. (41) and (39) we find the 2nd order approximation for the cross-sectionnally averaged concentration $C$:

[TABLE]

Eq. (40) is exactly the Taylor dispersion equation.

V.3 Taylor Dispersion with solute absorption conditions at the wall

Invariant manifold method with absorption in the steady uniform contracting channel

We consider the problem in the steady uniform case again but this time with absorption at the boundaries. The boundary condition at the boundaries of the channel with absorption writes:

[TABLE]

where $\gamma^{\prime}=\gamma\frac{a}{\kappa}$ and $\gamma$ is the surface absorption coefficient in $\mathrm{m.\,s^{-1}}$. In the following we drop the prime notation for $\gamma$ and keep it non-dimensional.

0th order expansion

As before we have $V^{0}=C$ the mean concentration over the cross-section.

1st order expansion

Now we solve for $V^{1}$. It verifies:

[TABLE]

We use the expression of $u(\eta)$, $V^{0}$:

[TABLE]

and integrating once and using the boundary conditions we find $G^{1}$:

[TABLE]

where we notice that $G^{1}$ now also contains some sink term due to absorption at the boundaries. Reporting in the integral expression for $V^{1}$, and using the fact that the cross-sectional average of $V^{1}$ has to vanish yields the following expression:

[TABLE]

2nd order expansion

Now we solve for $V^{2}$. It verifies:

[TABLE]

Now we integrate one first time and use the boundary condition on $V^{2}$ gives:

[TABLE]

To compare the expansion with other works it is interesting to search for the next order in the expansion. We report just the order 3 term here:

[TABLE]

And order 4 term:

[TABLE]

Heuristic derivation for absorption at the boundaries

From heuristic derivations of Taylor dispersion with absorption such as in Meigel and Alim (2018) we have222Note that in Meigel and Alim (2018) there is a typo in the formula given that we have updated here.

[TABLE]

This gives the coefficients plotted in the main text (gray) that are not as good to reproduce reality. This can be explained in the following way. As the heuristic expansion is a 2nd order expansion (and neglects all higher order terms), one should stop at 2nd order terms in $\gamma$ and in $\frac{\partial.}{\partial z}$, which gives:

[TABLE]

which is absolutely correct and corresponds to the invariant manifold expansion.

V.4 Taylor Dispersion in a straight channel with different flow profiles

V.4.1 Flow slip at the channel wall

We consider the same problem as before, in the simplest setting where $\frac{\partial a}{\partial t}=0$ and $\frac{\partial a}{\partial z}=0$. This time, we allow for slippage at the boundaries of the channel. In other words, instead of having the standard Poiseuille profile, we impose:

[TABLE]

where $b>0$ is a length called the slip length. If we solve for the flow in the channel we find then:

[TABLE]

with $\alpha=1+2b/a$, such that the velocity at the interface is non zero. We also have that the cross-sectionally averaged flow is:

[TABLE]

0th order expansion

The assumptions on the flow profile do not change the other boundary conditions such that $V^{0}=C$.

1st order expansion

Now we solve for $V^{1}$. It verifies:

[TABLE]

We use the expression of $u(\eta)$ and integrate once with respect to $\eta$:

[TABLE]

Note there that we did not add any integration constant. Automatically this constant is set to zero, otherwise, once divided by $\eta$, it would yield divergences in $\frac{\partial V^{1}}{\partial\eta}$. We now use the boundary condition in $\eta=1$, this yields the expression of $G^{1}$:

[TABLE]

We now use this to derive an expression for $V^{1}$. We find:

[TABLE]

and using the fact that the average over the cross section of $V^{1}$ is zero, $\int V^{1}\eta d\eta=0$ we find $B=-\frac{a^{2}}{12\kappa}U\alpha^{-1}\frac{\partial C}{\partial z}$, and thus:

[TABLE]

Step 2

Now we solve for $V^{2}$. It verifies:

[TABLE]

Now we use Eq. (38) and the expression for $u$; and integrate once:

Using the boundary condition on $V^{2}$ gives:

[TABLE]

[TABLE]

and with no further ado, assembling Eq. (48) and (50) yields the following second order approximate equation for the cross-sectionally averaged concentration $C$:

[TABLE]

Eq. (51) is exactly Eq. (14) of the main text.

V.4.2 Arbitrary flow profile

We consider the same problem as before, in the simplest setting where $\frac{\partial a}{\partial t}=0$ and $\frac{\partial a}{\partial z}=0$. The flow profile is some function of the coordinate $\eta$, $u(\eta)$. We do not impose any condition on it for now (except standard regularity). The orthogonal flow is still zero.

0th order expansion

This flow profile does not change the other boundary conditions such that $V^{0}=C$.

1st order expansion

Now we solve for $V^{1}$. It verifies:

[TABLE]

We use the expression of $u(\eta)$ and integrate once with respect to $\eta$:

[TABLE]

and using the boundary condition in $\eta=1$, this yields the expression of $G^{1}$:

[TABLE]

where $U=2\int_{0}^{1}u(x)xdx$ is the mean flow. We now use this to derive an expression for $V^{1}$. We find:

[TABLE]

At this stage let’s introduce $\tilde{u}(\eta)=\int_{0}^{\eta}\frac{1}{y}\int_{0}^{y}u(x)xdx$. We use the fact that the average over the cross section of $V^{1}$ is zero, $\int V^{1}\eta d\eta=0$ we find $B=\frac{a^{2}}{8\kappa}U\frac{\partial C}{\partial z}-\frac{a^{2}}{\kappa}\overline{\tilde{u}}\frac{\partial C}{\partial z}$, and thus:

[TABLE]

2nd order expansion

Now we solve for $V^{2}$. It verifies:

[TABLE]

Now we use Eq. (53) and integrate once:

[TABLE]

Using the boundary condition on $V^{2}$ after some straightforward calculus:

[TABLE]

where

[TABLE]

such that we get the following 2nd order approximate equation for the cross sectionnally averaged concentration $C$:

[TABLE]

Eq. (55) is exactly Eq. (15) in the main text.

V.4.3 Optimal flow profile

First, we make the change of variable: $\psi(\eta)=\frac{u(\eta)}{U}-1$. $\overline{\psi}=0$ and measures the deviation of the velocity profile $u$ to a plug flow. One can rewrite Eq. (54) as:

[TABLE]

There the space of integration is for all $x\geq(x^{\prime},x^{\prime\prime})$ such that we can also rewrite:

[TABLE]

We search for the optimal flow profile $v$, that maximizes $I[u]$. We will impose a few conditions. First, the average over the cross-section has to vanish:

[TABLE]

and second, we impose that the $L^{2}$ norm is bounded, namely that the deviation to the plug flow has a finite strength:

[TABLE]

where $\alpha=1/6$ is taken from the Poiseuille flow reference.

The problem then amounts to finding an extremal solution of the following lagrangian:

[TABLE]

To do so we search for $v$ such that for any function $h(x)$ we have $\mathcal{L}[\psi+h]-\mathcal{L}[\psi]=0$. This condition writes:

[TABLE]

This is equivalent to searching $\psi$ such that $\forall x$:

[TABLE]

and as this equation is true $\forall x$ it is true in particular for $x\neq 0$:

[TABLE]

and since $\int_{0}^{1}\psi(y)ydy=0$ we find:

[TABLE]

As this is true $\forall x$, it is true for the derivative of this function as well,

[TABLE]

multiplying by $x$ and performing the derivative one more time gives

[TABLE]

which is a second order partial differential equation. Assuming that $\psi(x)$ is sufficiently regular to be expandable as series: $\psi(x)=\sum_{n}a_{n}x^{n}$, one easily finds that the optimal solution $\psi$ is a Bessel function:

[TABLE]

To close the problem one has to find $\lambda_{1}$ such that Eq. (58) is verified and $a_{0}$ such that Eq. (59) is verified. $\lambda_{1}$ has multiple solutions:

[TABLE]

The first equality has multiple roots in $\lambda_{1}$. The largest root is $\lambda_{1}\simeq 0.068$ and actually corresponds to the optimal value, while the other roots corresponds to smaller values of $I[u]$. This value of $\lambda_{1}$ gives $a_{0}\simeq 1.43$ in the second equality.

V.5 Taylor Dispersion in a non uniform time contracting channel

We focus now on the case where the radius of the channel $a(z,t)$ depends on the axial position along the channel $z$ and time $t$. Before we can start solving the problem in this case, some ideas must be set up. We highlight at this point that the non dimensional variable $\eta(z,t)=\frac{r}{a(z,t)}$ depends on $z$ and time $t$.

Problem set up

We consider that the flow is quasi-Poiseuille, such that we have a similar expression for the longitudinal velocity $u$, but there is an axial component $v$ that can be computed easily by solving the zero divergence.

[TABLE]

Moreover, one must not forget that considering that the flow is incompressible yields the following useful relationship between the mean flow and the cross section area:

[TABLE]

The no solute flux boundary condition at the channel boundary still writes (at first order in $\frac{\partial a}{\partial z}$):

[TABLE]

Note that in this expression, $\frac{\partial a}{\partial t}$ does not appear because the condition no solute flux boundary condition may be written in the moving reference frame of the boundary.

0th and 1st order expansion

The 0th and 1st order of the expansion do not change as compared to the straight steady case and thus we may start off with: $V^{0}=C$, $G^{1}=-U\frac{\partial C}{\partial z}$ and $V^{1}=\frac{a^{2}}{4\kappa}U\frac{\partial C}{\partial z}\left[\eta^{2}-\frac{\eta^{4}}{2}-\frac{1}{3}\right]$.

2nd order expansion

Now we solve for $V^{2}$. It verifies:

[TABLE]

As compared to the previous case, we have to pay special attention to derivatives in space and time especially because the non-dimensional coordinate $\eta=r/a(z,t)$ is dependent in space and time as well:

[TABLE]

And also to $\frac{\partial V^{1}}{\partial t}$ that is non zero here:

[TABLE]

We now assemble everything and integrate once to find $G^{2}$ at the boundary. For simplicity, in the following derivation we include by a symbolic $\tilde{f}(\eta)=\left[\eta^{2}-\frac{\eta^{4}}{2}-\frac{1}{3}\right]$. All terms that are not multiplied by another function of $\eta$ are forgotten at the second line because they cancel out in the derivation:

[TABLE]

[TABLE]

Now we use the boundary condition:

[TABLE]

and thus:

[TABLE]

And thus if we assemble the parts to get the partial differential equation for $C$ we find:

[TABLE]

Note at this point that $\frac{\partial a^{2}U}{\partial z}=-2\frac{\partial a}{\partial t}a$, and thus one can rewrite the previous equation in the simpler form :

[TABLE]

V.6 Effect of pumping on dispersion

We start with the differential equation on $\rho$

[TABLE]

At second order in $\phi$ we have:

[TABLE]

with the phase velocity $v_{p}=\omega/k$ and

[TABLE]

Clearly the problem is periodic with $\tilde{U}$ and $\tilde{\kappa}$ periodic in $L=2\pi/k$. We are facing a one-dimensional advection-diffusion problem, for which the long time diffusion coefficient writes

[TABLE]

where $x(t)$ represents a stochastic variable, typically the position of the particle. The long time diffusion coefficient can be calculated either by the method of moments of first passage times Reimann et al. (2001) or by investigating probability distributions as in Guérin and Dean (2015). It is simply expressed as

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Since we work in the limit $\phi\ll 1$ (small deformations) it is possible to expand the exponentials present in $I_{+}$ and $I_{-}$. This allows to compute the integrals analytically and after some lengthy though easy computations, one finds:

[TABLE]

Exactly the result of the main paper.

The effective drift may also be obtained

[TABLE]

In this setting one finds at second order

[TABLE]

Note that this effective speed is taken in the moving referential, so one has to subtract the relative referential speed to get the drift in the original referential. Also, to get the appropriate effective drift in the original referential, one should add the mean input current from the boundary, which at order 2 writes $\frac{\phi^{2}}{2}\frac{\omega}{k}$ and therefore in the original referential one has

[TABLE]