Interactions enhance dispersion in fluctuating channels via emergent flows

Abstract

Understanding particle motion in narrow channels can guide progress in numerous applications, from filtration to vascular transport. Thermal or active fluctuations of fluid-filled channel walls can slow down or increase the dispersion of tracer particles via entropic trapping in the wall bulges or hydrodynamic flows induced by wall fluctuations, respectively. Previous studies concentrated primarily on the case of a single Brownian tracer. Here, we address what happens when there is a large ensemble of interacting Brownian tracers – a common situation in applications. Introducing repulsive interactions between tracer particles, while ignoring the presence of a background fluid, leads to an effective flow field. This flow field enhances tracer dispersion, a phenomenon reminiscent of that seen for single tracers in incompressible background fluid. We characterise the dispersion by the long-time diffusion coefficient of tracers numerically and analytically with a mean-field density functional analysis. We find a surprising effect where an increased particle density enhances the diffusion coefficient, challenging the notion that crowding effects tend to reduce diffusion. Here, inter-particle interactions push particles closer to the fluctuating channel walls. Interactions between the fluctuating wall and the now-nearby particles then drive particle mixing. Our mechanism is sufficiently general that we expect it to apply to various systems. In addition, our perturbation theory quantifies dispersion in generic advection–diffusion systems.

1. Introduction

Fluid transport in confined channels, and generally in porous structures, is relevant for a broad range of biological and industrial applications, from nutrient transport in vascular networks and microorganisms in soils (Bhattacharjee & Datta 2019; Tomkins, Hughes & Morris 2021), to improved filtration of liquids or gases, including modern desalination techniques and oil recovery (Werber, Osuji & Elimelech 2016; Marbach & Bocquet 2019). However, narrow channels present a number of challenging features to model, e.g. the predominance of surface effects, the importance of spatiotemporal fluctuations, as well as specific electrostatic response (Kavokine, Netz & Bocquet 2021). Significant progress in the last decade has improved our understanding of transport features in such porous environments, and we briefly review below the effects of (i) spatial and (ii) temporal fluctuations of the confining environment, as well as (iii) crowding effects due to interactions.

First, purely spatial corrugations of the confining channel reduce the long-time diffusion coefficient of isolated particles. In fact, random crossings of channel constrictions are rare events that impede overall transport: the constrictions form effective entropic barriers (Zwanzig 1992). This effect is well captured by the so-called Fick–Jacobs formalism (Jacobs 1935; Reguera & Rubí 2001; Kalinay & Percus 2006; Rubi 2019), which reduces the problem to an effectively one-dimensional one. The Fick–Jacobs formalism has been extended recently to arbitrary channel geometries (Dagdug, García-Chung & Chacón-Acosta 2016; Chávez, Chacón-Acosta & Dagdug 2018). Recent work suggests that the approach is also adapted to study fluid flow through biological membranes (Arango-Restrepo & Rubi 2020). Such entropic contributions induce significant corrections to transport in microfluidic (Yang et al. 2017) or biological (Rubí et al. 2017) channels.

Second, and of importance in several applications, confining channels are often not static but fluctuate in time, due to either thermal agitation or an external forcing. Molecular dynamics simulations found, early on, enhanced diffusion of gas in microporous materials if the thermal vibrations of the material are accounted for (Leroy, Rousseau & Fuchs 2004; Haldoupis et al. 2012), and more recently, enhanced water diffusion in solid state pores via phonon-fluid coupling (Ma et al. 2015, 2016; Cao, Wang & Ma 2019; Noh & Aluru 2021). Recent theoretical work suggests that enhanced diffusion is universally due to longitudinally fluctuating fluid flows, driven by fluctuations of the channel walls, that convect the tracer particles (Marbach, Dean & Bocquet 2018). The mechanism is thus reminiscent of Taylor–Aris dispersion, where the long-time diffusion coefficient is enhanced by a cross-sectionally inhomogeneous fluid flow profile (in the initially studied case of Poiseuille flow) (Taylor 1953; Aris 1956; Hydon & Pedley 1993). When the characteristic time scale of wall fluctuations is smaller (larger) than the time scale to diffuse across typical constrictions, diffusion is enhanced (decreased), a criterion that was verified numerically (Chakrabarti & Saintillan 2020) and in experiments in fluctuating porous matrices (Sarfati, Calderon & Schwartz 2021). Such enhancement of dispersion properties is relevant in a number of biological contexts, such as in blood vessels (Masri, Puelz & Riviere 2021), slime mould vasculature (Marbach & Alim 2019), the gut (Codutti, Cremer & Alim 2022), and near molecular motors (Evans, Krause & Feringa 2021), and are of general relevance in plants (Tomkins et al. 2021).

Third, in addition to the fluctuating confinement, tracer particles are often not isolated and interact with other particles or molecules that are also diffusing in the medium. The current paradigm is that such crowded environments, in an open domain, tend to slow down diffusion at equilibrium (Lekkerkerker & Dhont 1984; Lowen & Szamel 1993; Dean & Lefèvre 2004). However, subtle effects may arise if the crowded medium is driven out of equilibrium (Zaccone et al. 2011; Lappala, Zaccone & Terentjev 2013). For example, the diffusion coefficient of a tracer driven by an external force may be enhanced at low density and high forcing (Bénichou et al. 2013, 2018; Démery, Bénichou & Jacquin 2014; Illien et al. 2018). Typically, a trade-off is observed between, on the one hand, the increased diffusion induced by faster exploration of space thanks to the driving force, and on the other hand, the decreased diffusion due to spatial constrictions induced by the confining media (Illien et al. 2018). Such effects may be exploited for active microrheology within spatially corrugated channels (Puertas, Malgaretti & Pagonabarraga 2018; Malgaretti, Puertas & Pagonabarraga 2022).

Given its obvious importance, especially in biological systems, the coupling between the effects of fluctuating channels and inter-particle interactions has received surprisingly little attention. Since temporal channel fluctuations increase transport coefficients, and since inter-particle interactions, or crowding effects, generally decrease diffusion, it is natural to ask if one can predict which effect dominates. Even though recent numerical work showed that diffusion enhancement could be obtained with increased particle density in microporous matrices (Obliger et al. 2019), sometimes even exhibiting a maximum (Pireddu et al. 2019) at a certain density, it was also noticed that the effect can depend strongly on the precise details of the system understudy (Obliger et al. 2023). These results thus call for a general theoretical investigation.

In this paper, we investigate the motion of diffusing particles with repulsive interactions in a confined and fluctuating channel (see figure 1a), which is essentially a spatially periodic profile moving with constant velocity $v_{wall}$. We perform numerical simulations to quantify the effective long-time self-diffusion coefficient $D_{eff}$ and the effective long-time drift $V_{eff}$ of particles. We explain their behaviour for a broad range of interaction strengths between particles and fluctuating channel speeds $v_{wall} = \omega _0/k_0$, where $k_0$ and $\omega _0$ are the wavenumber and frequency of the wall fluctuations. Using perturbation theory, we obtain simple analytical predictions for ${V_{eff}}$ and $D_{eff}$ that are in excellent agreement with simulations.

Figure 1. Set-up to study the transport of tracers in an interacting system with fluctuating boundaries. (a) Tracer particles (yellow) perform a random walk within this wiggling environment, here represented by a top, moving wall (red arrows). (b) We consider pairwise, soft, repulsive interactions between particles, and we vary the interaction strength $\alpha$ and the particle number density $\rho _0$ to investigate more or less crowded environments.

The paper is organised as follows. To start, in § 2, we consider a simple ideal gas in the fluctuating channel. We find that the long-time diffusion coefficient $D_{eff}$ exhibits a maximum with respect to the wall phase velocity $v_{wall}$, corresponding to maximal enhancement of diffusion. This result should be contrasted with that of the diffusion of tracer particles in incompressible fluids that exhibits a monotonic increase with the wall phase velocity $v_{wall}$ (Marbach et al. 2018). The non-monotonic behaviour seen in the case studied here originates from the interplay between diffusion and advection due to the wall that increases long-time diffusion only when the diffusive time scale and the advection time scale are comparable. In § 3, we then consider soft-core interactions between particles (see figure 1b). We find that $D_{eff}$ can be enhanced further with increasing repulsive interactions. Here, repulsive interactions play a role in generating a more uniform distribution of particles in the channel, even in the vicinity of rapidly fluctuating bulges. Eventually, increased wall–particle collisions, caused by now-nearby particles, enhance dispersion. This behaviour is reminiscent of the mechanism of enhanced tracer diffusion in a fluctuating channel filled with an incompressible fluid (Marbach et al. 2018), and indeed in § 4, we show that remarkably, analytically and numerically, transport coefficients converge to a universal incompressible fluid regime in the limit of strong repulsive interactions. Finally, in § 5 we discuss how these mechanisms and techniques may be used further to investigate more complex situations of transport in fluctuating confined environments.

2. Transport of ideal (isolated) tracers in a fluctuating channel

2.1. Simulation results

We start by exploring the motion of tracers in a fluctuating channel, where the environment is not crowded, i.e. where tracer particles are sufficiently far away from each other that we can consider that they do not interact. The tracer particles perform a random walk with diffusion coefficient $D_0$. We refer to this case as the ‘ideal gas’ regime. Throughout this study, we will assume that the impact of the fluctuating channel boundaries on the velocity field of the supporting fluid is negligible. This is the case of particles embedded in a highly compressible fluid or gas, for which the mean free path is much smaller than any relevant length scale in the system. Examples of such compressible systems in biology include airflow in the pulsating alveoli in the lungs (Sznitman et al. 2007; Dong, Yang & Zhu 2022). Alternatively, this situation can also be achieved considering that the boundaries are not implemented by hard walls but rather by potential barriers, such as (fluctuating) electrostatic potentials for charged tracers. The impact of fluctuating boundaries on the fluid's velocity field has already been characterised in the Stokes flow regime (Marbach et al. 2018).

We simulate the motion of these non-interacting particles in a fluctuating channel, described through a fluctuating wall at height $h(x,t) = H + h_0 \cos ( k_0 x - \omega _0 t )$. The shape of the fluctuating interface is chosen to be sinusoidal, as a generic interface profile can be decomposed in terms of plane waves. The presence of the walls is incorporated via a soft boundary potential acting on the particles. We track the motion of particles over long times, and evaluate their effective long-time self-diffusion coefficient $D_{eff}$ and mean drift $V_{eff}$ along the main channel axis $x$ (see Appendix A for details).

Previous work on incompressible fluids has established that a relevant parameter to analyse the system is given by the Péclet number characterising the fluctuations,

(2.1)\begin{equation} Pe = \frac{\omega_0}{D_0 k_0^2} = \frac{\tau_{diff}}{\tau_{wall}}, \end{equation}

which compares the time scale to diffuse across the length of a channel corrugation $\tau _{diff} = 1/D_0 k_0^2$ to the period of the channel fluctuations $\tau _{wall} = 1/\omega _0$ (Marbach et al. 2018). In simulations, we therefore fix the typical channel corrugation length $L = 2 {\rm \pi}/ k_0$ and vary the wall phase velocity $v_{wall} = \omega _0/k_0$. All parameters are expressed in terms of a time unit $\tau _0$ and a length unit $\ell _0$, which are arbitrary. We present the results (yellow dots) in figures 2(a,b) for $D_{eff}$ and ${V_{eff}}$ with increasing $Pe$.

Figure 2. Transport properties of non-interacting tracers in a periodic driven channel; comparison between a compressible (ideal gas) and an incompressible supporting fluid. (a) Long-time longitudinal diffusion coefficient and (b) drift of particles, with $Pe = \omega _0/D_0 k_0^2 =v_{wall}L/(2{\rm \pi} D_0)$. In (a) and (b), error bars corresponding to one standard over 10 independent runs are smaller than point size, except for (b) small $Pe$, since $V_{eff} \simeq v_{wall} \sim 0.01 \ell _0/\tau _0$ is small compared to the noise level. Ideal gas theory corresponds to (2.18) and (2.19), and incompressible fluid to (2.24) and (2.25). Stationary density profiles of particles for $v_{wall}=0.1 \ell _0/\tau _0$ (or $Pe \simeq 3$) for (c) the ideal gas and (d) an incompressible fluid, for $v_{wall}=0.5 \ell _0/\tau _0$ ($Pe \simeq 16$). Blue arrows show the velocity field in the channel, represented in the referential of the moving wall, with length proportional to magnitude (arbitrary scale). The colour scale for the density profiles is shared in (c) and (d), and yellow (purple) regions indicate regions of high (low) density. Other numerical parameters are $L=200\ell _0$, $H=12 \ell _0$, $h_0=3\ell _0$ and $D_0=1 \ell _0^2/\tau _0$.

Interestingly, the long-time diffusion coefficient $D_{eff}$ exhibits a non-monotonic dependence on the Péclet number $Pe$ (yellow dots in figure 2a), which we interpret in the following way. At low Péclet numbers $Pe \ll 1$, particles move much faster than the wall, $\tau _{wall} \gg \tau _{diff}$, hence the wall appears frozen. Particles therefore spend time exploring the bulges before escaping the constrictions, and their diffusion coefficient is consequently decreased, $D_{eff} \leq D_0$. This effect is well captured by the Fick–Jacobs approximation, where transport is reduced due to constrictions acting as effective entropy barriers (Jacobs 1935; Zwanzig 1992; Reguera & Rubí 2001; Kalinay & Percus 2006; Burada et al. 2007; Mangeat, Guérin & Dean 2017; Marbach et al. 2018; Rubi 2019). At intermediate Péclet numbers $Pe \gtrsim 1$, $\tau _{wall} \simeq \tau _{diff}$ and the moving boundaries enhance particle motion: particles collide with the moving corrugated walls, which increases their diffusion coefficient overall, $D_{eff} \geq D_0$. At high Péclet numbers $Pe \gtrsim 10$, $\tau _{wall} \ll \tau _{diff}$ and the wall moves so fast that particles no longer have time to enter the bulges; they therefore behave as though they were not seeing any corrugation, the effective diffusion coefficient is unchanged and equals the bare diffusion, $D_{eff} = D_0$. In other words, from the point of view of the particles, the channel is effectively flat with height given by its minimum height.

To further understand the behaviour at high Péclet numbers, we plot the average density distribution $\rho (x,y)$ of particles within the channel (the average being over noise), in the reference frame where the channel profile is stationary, for a large value of $v_{wall}$ (figure 2c). We observe that particles accumulate at the constriction. Indeed, in this reference frame, the channel constriction acts as a bottleneck for transport. This can be seen as an inverse Bernoulli effect: in the channel's frame of reference, the particle flux $\rho (x) h(x) v_{wall}$ is necessarily constant, where $\rho (x)$ is the cross-sectional averaged density of particles. Hence the density (and not the velocity, as in the Bernoulli effect) is largest where the channel is constricted, $\rho (x) \propto 1/h(x)$. Eventually, since the particles explore only low vertical coordinates $y \lesssim H - h_0$, they no longer collide with the moving wall, hence their diffusion coefficient is unchanged.

The effective particle drift ${V_{eff}}$ decreases monotonically with the Péclet number (yellow dots in figure 2(b), noting that the vertical axis is the normalised drift relative to the wall phase velocity, ${V_{eff}} /v_{wall}$). At small Péclet numbers, the density of particles is uniform in the channel, therefore the particles within the bulge are carried along in the same direction as the wall phase velocity. At higher Péclet numbers, the fraction of particles within the bulge decreases, as seen in figure 2(c). Since all particles accumulate within the bottleneck, they are no longer pushed by the moving wall, hence ${V_{eff}} /v_{wall} \rightarrow 0$.

2.2. Analytic theory to account for transport in fluctuating channels

To account for this broad behaviour, we build a general analytic theory that reproduces these effects. With the goal of making a pedagogical introduction to our perturbation method, which closely follows that which was described only briefly in Marbach et al. (2018), we devote this subsection to a detailed explanation.

2.2.1. Constitutive equations

Brownian tracer particles evolve in a fluctuating environment, confined in the $y$ direction between $y = 0$ and $y = h(x,t)$, and infinitely long in the $x$ direction. Compared to the simulation, we here make no assumption on the form of $h(x,t)$, which may describe thermal motion of the wall (determined e.g. through a Hamiltonian characterising the flexibility of the interface) or active motion (driven by an external source of energy, as in our simulations where fluctuations are imposed). The probability density function $\rho _{tr}(x,y,t)$ to find a tracer particle at position $(x,y)$ at time $t$ obeys the Fokker–Planck equation

(2.2)\begin{equation} \frac{\partial{\rho_{tr}(x,y,t)}}{\partial t}+ \boldsymbol{ \nabla} \boldsymbol{\cdot} \boldsymbol{j} (x,y,t) = 0, \end{equation}

with $\boldsymbol {j}(x,y,t) \equiv - D_0\,\boldsymbol { \nabla } \rho _{tr}(x,y,t) + \boldsymbol {u} (x,y,t)\,\rho _{tr}(x,y,t)$, and where we have assumed that the tracer's diffusion coefficient $D_0$ is uniform in space and that the tracer is advected by the field $\boldsymbol {u} (x,y,t)$, which can have both potential and non-potential components. In the case of ideal non-interacting particles, the underlying supporting fluid does not move due to the moving interface, and there are no interactions, $\boldsymbol {u} (x,y,t) \equiv 0$. For now, we keep $\boldsymbol {u}(x,y,t)$ arbitrary as this will be useful in the incompressible and interacting regimes.

2.2.2. Boundary conditions

We impose no flux boundary conditions at both surfaces. This means that the projection of the current on the direction normal to the lower surface boundary is zero, and that similarly, the projection of the current on the direction normal to the upper boundary in the frame moving with speed $\boldsymbol {U}_\mathrm {wall}$, where the surface is static, is zero. The boundary conditions thus read

(2.3)\begin{gather} j_y(x,0,t) = 0, \end{gather}

(2.4)\begin{gather}\boldsymbol{j} (x,h(x,t),t) \boldsymbol{\cdot} \boldsymbol{n} =- \rho (x,h(x,t),t)\,\boldsymbol{U}_{wall} \boldsymbol{\cdot} \boldsymbol{n}. \end{gather}

Here, $\boldsymbol {n}$ is the outward normal to the interface, and $\boldsymbol {U}_{wall} = (0, \partial _t h)$ is the velocity of the wall, considering that the wall atoms move vertically about their average position, similarly to phonon modes that propagate on material structures (Chaikin, Lubensky & Witten 1995) or peristaltic motion of vasculature (Alim et al. 2013; Marbach & Alim 2019). We derive these boundary conditions from first principles in § 1 of the supplementary material available at https://doi.org/10.1017/jfm.2023.640.

2.2.3. Simplified longitudinal equation for the probability distribution

To make progress on the long-time behaviour of (2.2), we place ourselves within the lubrication approximation. This means that we consider the typical corrugation length $L$ to be much bigger than the average channel height, $\langle h(x,t) \rangle = H$, $H \ll L$, itself much bigger than the amplitude of the channel fluctuations, i.e. $\sqrt {\langle (h - H)^2 \rangle }\propto h_0 \ll H$. Within the lubrication approximation, the outward normal to the interface in (2.4) simplifies to $\boldsymbol {n} \simeq (\partial _x h,1)$. In this framework, the particle density relaxes much faster on the vertical direction $y$ than on the longitudinal direction $x$. At low Péclet number, this should notably yield a vertically uniform particle density. Figure 2(c) shows that the density profile is indeed independent of $y$ in our simulations. Also, in Appendix B, we show that even though the lubrication approximation should fail for e.g. larger channel heights $H$, the results derived below remain surprisingly robust.

We therefore look for an approximate evolution equation on the probability distribution integrated vertically, or marginal distribution, $p_{tr}(x,t) = \int _0^{h(x,t)} \rho _{tr}(x,y,t) \, \text {d} y$. Taking the time derivative of $p_{tr}(x,t)$ and using (2.2) yields

(2.5) \begin{align} \frac{\partial{p_{tr}(x,t)}}{\partial t} &= \frac{\partial{h(x,t)}}{\partial t}\,\rho_{tr}(x,h(x,t),t)- \int_0^{h(x,t)} \left( \frac{\partial{j_y(x,y,t)}}{\partial y} + \frac{\partial{j_x(x,y,t)}}{\partial x} \right) \text{d} y \nonumber\\ &= \frac{\partial{h(x,t)}}{\partial t}\,\rho_{tr}(x,h(x,t),t) - j_y(x,h(x,t),t) + j_y(x,0,t) \nonumber\\ &\quad + \frac{\partial{h(x,t)}}{\partial x}\,j_x(x,h(x,t),t) - \frac{\partial{}}{\partial x} \left(\int_0^{h(x,t)} j_x(x,y,t) \,\text{d} y\right), \end{align}

where we used the simple and useful relation

(2.6)\begin{equation} \frac{\partial{}}{\partial x} \left( \int_0^{h(x,t)} f(x,y,t) \,\text{d} y\right) = \frac{\partial{h(x,t)}}{\partial x}\,f(x,h(x,t)) + \int_0^{h(x,t)} \frac{\partial{f(x,y,t)}}{\partial x} \,\text{d} y \end{equation}

for any function $f$. Equation (2.5) can be simplified remarkably by using the no-flux boundary conditions (2.4) and (2.3), which lead to all the surface terms vanishing, and we find

(2.7)\begin{equation} \frac{\partial{p_{tr}(x,t)}}{\partial t} =- \frac{\partial{}}{\partial x} \left(\int_0^{h(x,t)} j_x(x,y,t) \,\text{d} y\right). \end{equation}

We look for a closed equation on $p_{tr}(x,t)$. Writing

(2.8)\begin{equation} j_x(x,y,t) =- D_0\,\frac{\partial{\rho_{tr}(x,y,t)}}{\partial x} + u_x\, \rho_{tr}(x,y,t) \end{equation}

explicitly, and using (2.6) again, we find

(2.9)\begin{align} \frac{\partial{p_{tr}(x,t)}}{\partial t}& =- \frac{\partial{}}{\partial x} \left(- D_0\,\frac{\partial p_{tr}(x,t)}{\partial x} + D_0\,\frac{\partial{h(x,t)}}{\partial x}\, \rho_{tr}(x,h(x,t),t)\right.\nonumber\\ &\quad \left.{}+ \int_0^{h(x,t)}u_x\,\rho_{tr}(x,y,t) \,\text{d} y \right), \end{align}

where the last terms still do not depend explicitly on the marginal $p_{tr}$. We will therefore make the common first-order approximation that since the probability distribution profile is nearly uniform vertically, we may assume $\rho _{tr}(x,y,t) \sim p_{tr}(x,t)/h(x,t)$. We then obtain the following Fokker–Planck equation on the marginal distribution, at lowest order in $\rho _{tr}(x,y,t) - p_{tr}(x,t)/h(x,t)$:

(2.10)\begin{equation} \frac{\partial{p_{tr}}}{\partial t} =- \frac{\partial{}}{\partial x} \left(-D_0\,\frac{\partial{p_{tr}}}{\partial x} + D_0 p_{tr}\,\frac{\partial{\ln h}}{\partial x} + \overline{u_x} p_{tr} \right), \end{equation}

where $\overline {u_x} = ({1}/{h}) \int _0^{h(x,t)} u_x \,\text {d} y$ is the vertically averaged longitudinal drift. This final step can be made more rigorous using a centre manifold expansion; see Mercer & Roberts (1990, 1994) and Marbach & Alim (2019). Equation (2.10) can alternatively be obtained by starting with the ansatz $\rho _{tr}(x,y,t) \sim p_{tr}(x,t)/h(x,t)$ at the level of (2.2), and making use of the boundary conditions. However, the derivation that we present here has the advantage of being approximation-free until (2.9).

Equation (2.10) clearly simplifies the initial problem, and it is sufficient to study the long-time behaviour of $p_{tr}$ to obtain the long-time diffusion coefficient $D_{eff}$ and drift ${V_{eff}}$. When $\overline {u_x} \equiv 0$, (2.10) corresponds exactly with the Fick–Jacobs equation (Jacobs 1935; Zwanzig 1992; Reguera & Rubí 2001), which describes the motion of (non-interacting) particles in spatially varying but time independent channels $y \equiv h(x)$. Interestingly, the Fick–Jacobs equation is therefore also valid for fluctuating channels in time, $h(x,t)$, regardless of the functional shape of $h(x,t)$, the only assumption being the lubrication approximation. We note that (2.10) is also consistent with the case of a background incompressible fluid, when $\overline {u_x}$ is the cross-sectionally averaged fluid velocity (Eq. (1) of Marbach et al. 2018).

2.2.4. Perturbation theory to obtain the long-time transport behaviour

Our goal is now to obtain the long-time transport behaviour of particles from the simplified, longitudinal evolution equation (2.10). We seek an approximate long-time equation for the marginal distribution $p_{eff}(x,t) = p_{tr}(x, t \rightarrow \infty )$ of the form

(2.11)\begin{equation} \frac{\partial{p_{eff}(x,t)}}{\partial t} = D_{eff}\,\frac{\partial^2 p_{eff}(x,t)}{\partial x^2} - V_{eff}\,\frac{\partial{p_{eff}(x,t)}}{\partial x} + \delta(x)\,\delta(t), \end{equation}

such that we can naturally read off the long-time diffusion $D_{eff}$ and drift $V_{eff}$, and where we assumed, without loss of generality, that the particle was initially at the centre of the domain.

To keep the calculations general, we rewrite (2.10) as a Fokker–Planck equation

(2.12)\begin{equation} \frac{\partial{p_{tr}(x,t)}}{\partial t} =- \frac{\partial{}}{\partial x} \left( - D_0\,\frac{\partial{p_{tr}(x,t)}}{\partial x} + v(x,t)\,p_{tr}(x,t) \right) + \delta(x)\,\delta(t), \end{equation}

where $v(x,t)$ is a general advection coefficient (for the ideal gas case, $v = D_0 \partial _x h / h$). Here, we consider that $v(x,t) = O(\varepsilon )$ is a fluctuating perturbation. Its average over realisations of the noise vanishes, $\langle v(x,t)\rangle = 0$, and its fluctuations are described in Fourier space by a spectrum $S(k,\omega )$ as

(2.13)\begin{equation} \langle \tilde{v}(k,\omega)\,\tilde{v}(k',\omega') \rangle = S(k,\omega)\,(2{\rm \pi})^2\,\delta(k + k')\,\delta(\omega + \omega'), \end{equation}

where $k$ and $\omega$ are the wavenumber and frequency, respectively, and the Fourier transform $\tilde {v}(k,\omega )$ of $v(x,t)$ is defined in (C1). Performing a perturbation development to solve (2.12) on the small parameter $\varepsilon$, we find (see Appendix C)

(2.14)\begin{gather} D_{eff} = D_0 \left(1- \varepsilon^2 \int \frac{\text{d} k \,\text{d} \omega}{(2{\rm \pi})^2}\,\frac{k^2( D_0^2 k^4 - 3 \omega^2)}{( D_0^2 k^4+ \omega^2 )^2}\,S(k,\omega)\right), \end{gather}

(2.15)\begin{gather}V_{eff} = \varepsilon^2 \int \frac{\text{d} k \,\text{d} \omega}{(2{\rm \pi}^2)}\,\frac{\omega}{k} \frac{k^2}{D_0^2k^4 + \omega^2}\,S(k,\omega). \end{gather}

Equations (2.14) and (2.15) are two of the main results of this work. They predict the long-time transport properties of particles within a fluctuating channel with arbitrary fluctuating local drift $v(x,t)$. The results can be applied regardless of the nature of the fluctuations, be they thermal or non-equilibrium, and regardless of the strength of the interactions between particles, and between particles and the supporting fluid. In essence, we generalise the results of Marbach et al. (2018), which were valid only for an incompressible supporting fluid.

2.3. Applications of the theory to the periodic channel

2.3.1. Ideal gas

To use the above analytic framework to describe our simulations, we now take the case of the periodic travelling wave $h(x,t) = H + h_0 \cos ( 2{\rm \pi} (x - v_{wall} t ) /L ) = H +$ $h_0 \cos ( k_0 x - \omega _0 t)$ on the surface. The Fourier transform of $h - H$ is $\tilde {h}(k,\omega ) = (h_0/2)$ $(\delta (k+k_0)\,\delta (\omega - \omega _0) + \delta (k-k_0)\,\delta (\omega + \omega _0) )$. The local drift $v$ in the ideal gas case is at lowest order in $\varepsilon = h_0/H$,

(2.16)\begin{equation} v(x,t) = D_0\,\frac{1}{H}\,\frac{\partial{h}}{\partial x},\quad \text{and in Fourier space}\quad \tilde{v}(k,\omega) = {\rm i} D_0 k\,\frac{\tilde{h(k,\omega)}}{H}. \end{equation}

This yields a spectrum for the local drift as

(2.17)\begin{equation} S(k,\omega) = \frac{D_0^2 k^2}{(2{\rm \pi})^2}\,\frac{h_0^2}{4 H^2}\,( \delta(k+k_0)\,\delta(\omega - \omega_0) + \delta(k-k_0)\,\delta(\omega + \omega_0) ). \end{equation}

Plugging the expression for the spectrum in (2.14) and (2.15) yields the long-time transport coefficients in the ideal gas case:

(2.18)\begin{gather} \frac{D^{\textit{ideal gas}}_{eff}}{D_0} = 1 + \frac{h_0^2}{2 H^2}\,\frac{3Pe^2-1}{(Pe^2 + 1)^2}, \end{gather}

(2.19)\begin{gather}\frac{V^{\textit{ideal gas}}_{eff}}{v_{wall}} = \frac{h_0^2}{2 H^2}\,\frac{1}{1 + Pe^2}, \end{gather}

where we recall the expression of the Péclet number $Pe = \omega _0/D_0 k_0^2$.

We present plots of (2.18) and (2.19) as lines in figures 2(a,b), and find excellent agreement with simulations. This agreement is robust over a wide range of physical parameters – see figure 9. Analytically, we recover for $Pe \rightarrow 0$ (corresponding to fixed channels, $v_{wall} = 0$) the well-known entropic trapping result where $D_{eff} = D_0 ( 1- h_0^2/2H^2)$ (Zwanzig 1992; Jacobs 1935; Reguera & Rubí 2001; Marbach et al. 2018; Rubi 2019). The analytic computation recovers the non-monotonic behaviour of the diffusion coefficient for intermediate $Pe$, and confirms that $D_{eff} = D_0$ at high $Pe\to \infty$. Finally, the amplitude of the correction for both the diffusion and the drift scales as $h_0^2/H^2$, which can be interpreted naturally and confirms the collision mechanism in bulges; indeed, the strength of the fluctuations scales as $h_0/H$ but only a fraction $h_0/H$ of particles lie within the bulge and take part in the enhancement of the diffusion coefficient or in the mean drift.

2.3.2. Comparison with transport in incompressible fluid

We now relate our results for the ideal gas to transport within incompressible fluids (Marbach et al. 2018). In that case, when the channel walls fluctuate, because the supporting fluid is incompressible, they induce fluctuations in the fluid's velocity field. Particles are thus advected by fluid flow. The channel walls here are perfectly slipping walls, which corresponds to the smooth boundary conditions considered for the gas particles, which have no specific lateral friction with the walls. The flow field $(u_x,u_y)$ can be derived in the low-Reynolds-number limit, which is the relevant limit to consider since we envision applications in microfluidics to nanofluidics. We find

(2.20)\begin{gather} u_x(x,y,t) = U_0(x,t), \end{gather}

(2.21)\begin{gather}u_y(x,y,t) = \frac{y}{h} \left( \frac{\partial{h(x,t)}}{\partial t} + U_0\,\frac{\partial{h(x,t)}}{\partial x} \right), \end{gather}

where $U_0(x,t)$ is the average fluid flow. We calculate $U_0(x,t)$ assuming peristaltic flow, i.e. that the flow is driven purely by channel fluctuations and that there is no mean pressure-driven flow (Marbach & Alim 2019; Chakrabarti & Saintillan 2020). The average pressure force on the fluid has to vanish, and we find ${U_0(x,t) \simeq v_{wall} (h_0/H)}$ ${\cos ( 2{\rm \pi} (x - v_{wall} t )/L )}$ at lowest order in $h_0/H$ (see § 2 of the supplementary material). The flow field is presented in figure 2(d) as blue arrows. As the channel moves towards the right-hand side, fluid mass in the right-hand side bulge has to swell out, consistently with the flow lines. Similarly, the bulge on the left-hand side opens up, allowing fluid flow to come in.

The advection term now has two contributions, one coming from the spatial inhomogeneities, and one coming from advection by fluid flow

(2.22)\begin{equation} v(x,t) = D_0\,\frac{1}{H}\,\frac{\partial{h(x,t)}}{\partial x} + U_0(x,t),\end{equation}

or $\tilde {v}(k,\omega ) = ( {\rm i} D_0 k + \omega /k )\,\tilde {h}(k,\omega )/H$ in Fourier space. The spectrum of the fluctuating drift is thus

(2.23)\begin{equation} S(k,\omega) = \frac{ D_0^2 k^2 + \omega^2/k^2 }{(2{\rm \pi})^2}\,\frac{h_0^2}{4 H^2}\,( \delta(k+k_0)\,\delta(\omega - \omega_0) + \delta(k-k_0)\,\delta(\omega + \omega_0) ). \end{equation}

We then obtain

(2.24)\begin{gather} \frac{D_{eff}^{\textit{inc. fluid}}}{D_0} = 1 + \frac{h_0^2}{2H^2}\,\frac{3 Pe^2 -1}{Pe^2 + 1}, \end{gather}

(2.25)\begin{gather}\frac{V_{eff}^{\textit{inc. fluid}}}{v_{wall}} = \frac{h_0^2}{2H^2}. \end{gather}

We perform numerical simulations where particles are advected by the flow field defined by (2.20) and (2.21). We present the numerical results as blue dots and the analytical results as blue lines in figures 2(a,b). We find perfect agreement between simulation and theory, confirming the analytical approach of Marbach et al. (2018).

In contrast with the ideal gas, when particles are surrounded by an incompressible fluid, $D_{eff}$ increases monotonically until it reaches a plateau at large $Pe$, and ${V_{eff}}$ stays constant regardless of $Pe$. In fact, at any value of $Pe$, and especially at high $Pe$, the density distribution of particles within the channel is uniform, as can be seen in figure 2(d). Therefore, even at high $Pe$, the population of particles lying within the bulges is pushed by the moving boundaries and increases both $D_{eff}$ and ${V_{eff}}$.

In this case it is natural to ask why the density distribution remains homogeneous along the channel. Looking closely at the velocity field (figure 2d), we see that particles are carried away from the bottleneck by the flow field, and into the bulges. The flow field therefore facilitates recirculation of accumulated particles. This naturally raises the question of how the supporting fluid's compressibility changes the transport properties of particles within fluctuating channels.

3. Interactions increase diffusion and drift in fluctuating channels

To investigate the impact of the compressibility of the supporting fluid, we introduce interactions between the particles (see figure 1b), and tune the interaction strength to vary the effective compressibility of the system.

3.1. Pairwise interactions and compressibility

We simulate the dynamics of interacting particles within a simple periodic fluctuating channel $h(x,t) = H + h_0\cos (k_0 x - \omega _0t)$, as in § 2. We use a pairwise interaction potential between particles,

(3.1)\begin{equation} \mathcal{V}_{int}(r) =\begin{cases} \alpha (r - d_c)^2 & \text{if } r < d_c,\\ 0 & \text{if } r \geq d_c, \end{cases} \end{equation}

where $r$ is the inter-particle distance, $d_c$ is a critical distance characterising the radius of the interaction (see figure 1b), and $\alpha$ is a spring constant that characterises the strength of the interaction. Note that $\alpha >0$ corresponds to repulsive interactions, while $\alpha <0$ corresponds to attractive interactions. We use a soft-core potential (instead of a hard-core potential as in e.g. Bénichou et al. 2013; Suárez, Hoyuelos & Mártin 2015) as it facilitates numerical integration over long time scales, which is necessary to obtain reliable statistics on the diffusion coefficient. Soft-core potentials are used commonly and also simplify analytic treatments (Pàmies, Cacciuto & Frenkel 2009; Démery et al. 2014; Antonov, Ryabov & Maass 2021). We will show later that the numerical results are well reproduced by our theory, whose predictions are robust to strong changes of the interaction potential. The mix of interacting particles and the surrounding (compressible) fluid forms a fluid of interacting particles, and we study its properties below.

The choice of interactions is well adapted to tune the compressibility $\chi ^{\star }_T$ of the fluid of interacting particles, and hence probe different compressibility regimes, in between the ideal gas and incompressible fluid cases. Indeed, $\chi ^{\star }_T$ is related to the structure factor $S_f(k)$ of the gas at zero wavelength (Hansen & McDonald 2013):

(3.2)\begin{equation} \chi^{{\star}}_{{T}}(\rho_0,\alpha) = \chi_T^{id}\lim_{k \rightarrow 0} S_f(k), \end{equation}

where $\chi _T^{id} = 1/(\rho _0 k_BT)$ is the compressibility of the ideal gas, which is infinite when $\rho _0 \rightarrow 0$. In the so-called and broadly used random phase approximation, one can relate the structure factor $S_f(k)$ to the interaction potential as

(3.3)\begin{equation} S_f(k) = \frac{1}{1 - \rho_0 \left( - \dfrac{1}{k_BT}\,\tilde{\mathcal{V}}_{int}(k)\right)}. \end{equation}

We calculate

(3.4)\begin{equation} \lim_{k\rightarrow 0} \tilde{\mathcal{V}}_{int}(k) = \iint \mathcal{V}_{int}(x,y) \,\text{d} x\, \text{d} y = \frac{\rm \pi}{6}\,\alpha d_c^4. \end{equation}

We therefore obtain the compressibility of the fluid of interacting particles as

(3.5)\begin{equation} \chi^{{\star}}_T(\rho_0,\alpha) =\frac{ \chi_T^{id}}{1 + \dfrac{\rm \pi}{6}\,\dfrac{\rho_0 \alpha d_c^4}{k_B T}} \equiv \frac{ \chi_T^{id}}{1 + \dfrac{E_0}{k_B T}}, \end{equation}

where we introduced $E_0 = ({\rm \pi} /6) \alpha d_c^4\rho _0$, which can be understood as the energy contribution of interactions contained in a typical volume $d_c^d$, with $d$ the dimension of space (here, $d=2$). For $\alpha <0$, one has $E_0<0$, and the fluid of particles ends up with a higher compressibility than the one of the ideal gas. In what follows, we will focus mainly on repulsive interactions ($\alpha >0$), but we keep in mind that our analytic computation does not assume the sign of $\alpha$. For $\alpha >0$, the compressibility $\chi ^{\star }_T(\rho _0,\alpha )$ decreases with the gas density $\rho _0$ and with increasing interaction strength. Therefore, either of the parameters $\rho _0$ and $\alpha$ is a good candidate to probe the intermediate regime between the ideal gas and the incompressible fluid for which $\chi _T^\star \to 0$. In the next subsection, we explore varying values of the density $\rho _0$, and we will find similar results when varying $\alpha$ (reported in Appendix D).

Finally, it is known that the long-time self-diffusion coefficient of interacting particles, in the bulk, i.e. in an open domain, $\overline{D}_0(\rho _0)$, is decreased in a crowded environment in equilibrium, compared to the infinitely dilute case (Démery et al. 2014) (although it may be non-monotonic at high densities as the energy landscape is flattened when the local density is large for soft interactions; see figure 7). We therefore first perform simulations within fixed flat walls ($h_0 = 0$, $v_{wall} = 0$), where the system is in equilibrium at temperature $T$. The confinement plays no particular role in that case because our channels are wide enough compared to the typical size of a particle ($d_c \ll H$). We measure the diffusion coefficients $\overline{D}_0(\rho _0)$ for various particle densities in the channel $\rho _0$ and interaction strengths $\alpha$, and our results agree well with existing theories (Démery et al. 2014; see Appendix A.4). We can thus now measure changes in the self-diffusion coefficient when there are fluctuations relative to this bulk value $\overline{D}_0(\rho _0)$.

3.2. Increasing interactions enhances diffusion and drift

We perform simulations of interacting particles in fluctuating channels ($h_0 = H/4$, varying $v_{wall} > 0$). We present our numerical results for the long-time self-diffusion coefficient $D_{eff}/\overline{D}_0(\rho _0)$ and mean drift ${V_{eff}} /v_{wall}$ with increasing particle density $\rho _0$ in figures 3(a,b). We find striking variations with increasing density. At low $Pe$, systems of interacting particles differ very little from systems with no interactions: we still observe entropic slow-down. Interestingly, at intermediate $Pe$, the enhancement of the diffusion coefficient increases with particle density $\rho _0$. The turnaround point at a critical value of $Pe$ increases also with increasing particle density. Similarly, mean particle drift is significant at larger values of $Pe$, the more so with increasing particle density. Increasing particle density thus does appear to bridge the gap between the ideal gas case and the incompressible fluid.

Figure 3. Transport properties of a fluid of interacting particles: (a) effective diffusion and (b) mean drift, for a compressible fluid of soft-core interacting Brownian particles with varying density $\rho _0$. Numerical parameters used here are similar to those in figure 2, with additionally $L=200 d_c$, $d_c = 2^{1/6}\ell _0 \simeq 1.12\ell _0$ and $\alpha =1 k_B T/\ell _0^2$. Error bars correspond to one standard deviation over 10 independent runs. Error bars for lower-density values are larger due to a smaller number of tracked particles (see table 1). Theory curves for the ideal gas and incompressible fluid are the same as for figures 2(a,b). Theory curves for the fluid of interacting particles correspond to (4.19) and (4.20).

3.3. Increasing interactions impacts the density profile

To understand the behaviour of $D_{eff}$ and ${V_{eff}}$, we need to understand first how particles rearrange due to inter-particle forces. Typically, a tracer particle will be forced down density gradients because of repulsive pairwise interactions. We therefore measure and analyse the particle density distributions in the channel for various average densities (see figures 4a,b). First, we find that the particle distribution is uniform in the vertical direction, as expected within the lubrication approximation: particles diffuse sufficiently fast on the vertical axis compared to all other relevant time scales. Second, compared to the ideal gas case at the same wall speed (figure 2(c) versus figure 4(a)), the particle distribution is also quite homogeneous on the horizontal direction, much like in the incompressible case (figure 2d). However, at higher wall velocities (figure 4b), particles accumulate again at the bottleneck region. Such accumulation in narrow channels is also observed in simulations of driven, interacting tracers in corrugated channels (Suárez et al. 2015; Suárez, Hoyuelos & Mártin 2016), which share a similar geometry in the reference frame where the wall is static.

Figure 4. Particle distribution within the channel for varying particle densities. (a,b) Particle distribution in a fluid of interacting particles at high density $\rho _0=10 \ell _0^{-2}$: (a) for intermediate wall velocity $v_{wall}=0.1 \ell _0/\tau _0$, and (b) for high wall velocity $v_{wall}=0.5 \ell _0/\tau _0$. The colour scale for the density profiles is shared between (a) and (b), and yellow (purple) regions indicate regions of high (low) density. (c,d) Marginal density $p(x)=\int \rho (x,y) \,\text {d} y$ of particles in the channel, calculated by integrating the two-dimensional density profile in (a) over vertical slabs as with the dashed grey box: (c) $p(x)$ for different values of $v_{wall}$, and (d) for different values of particle density $\rho _0$. Numerical parameters used here are similar to those in figure 3. Theory curves correspond to (4.10).

To further quantify the particle distributions, we study the marginal distribution profiles $p (x)=\int _{-\infty }^{\infty }\rho (x,y) \,\text {d} y$ along the channel, obtained in the simulations. As expected, $p(x)$ presents a peak, which indicates particles accumulating near the bottleneck. We find again that the profile $p(x)$ is more peaked with increasing wall velocity (figure 4c). It flattens out with increasing particle density $\rho _0$ (see figure 4d), indicating that density profiles indeed converge to the incompressible fluid case.

4. Analytic theory with pairwise interactions

4.1. Model for particle density profiles along the channel

To understand how the average density profile depends on interaction parameters, we expand our analytical theory. Our starting point is the equation governing the evolution of the density of particles $\rho (x,y,t)$, which can be written formally as (Dean 1996; Kawasaki 1998)

(4.1)\begin{equation} \frac{\partial{\rho}}{\partial t} = D_0\,\boldsymbol{ \nabla}^2 \rho + \frac{D_0}{k_B T}\,\boldsymbol{ \nabla} \boldsymbol{\cdot} ( \rho\,\boldsymbol{ \nabla} ( \mathcal{V}_{int} * \rho ) ) + \boldsymbol{\nabla} \boldsymbol{\cdot} ( \sqrt{2 D_0 \rho}\,\boldsymbol{ \xi} ),\end{equation}

where $\boldsymbol { \xi }$ is a two-component Gaussian white noise field, such that $\langle \xi _\mu (x,y,t) \rangle = 0$ and $\langle \xi _\mu (x,y,t)\,\xi _\nu (x',y',t') \rangle = \delta _{\mu \nu }\, \delta (x-x')\,\delta (y-y')\,\delta (t-t')$. Here, the convolution product is $(\mathcal {V}_{int} * \rho )(x,y) = \iint \text {d} x' \,\text {d} y'\,\mathcal {V}_{int}(\sqrt {x'^2 + y'^2})\,\rho (x- x', y-y')$, where $\mathcal {V}_{int}$ is any pairwise interaction potential. The term $D_0$ represents the microscopic diffusion constant of the individual Brownian particles. In an infinite domain and in the absence of interactions, the long-time diffusion constant of a tracer will be identical to $D_0$. The effects of interactions and geometry in the problem at hand here both play a role in modifying the late-time diffusion constant with respect to the microscopic one.

To make progress, we first decompose the density into average and fluctuating components:

(4.2)\begin{equation} \rho(x,y,t) = \bar{\rho}(x,y,t)+ \psi(x,y,t), \end{equation}

where $\psi (x,y,t)$ is the noise field such that $\langle \psi (x,y,t) \rangle = 0$. The noise perturbation $\psi (x,y,t)$ can be inferred, as in Démery et al. (2014), by assuming a large, flat, fixed environment and in the case where the background concentration $\bar {\rho }(x,y,t) = \rho _0$ is uniform. The short distance fluctuations in the density field can then be absorbed into a renormalization of the diffusion constant, which now depends on the average density. We will denote by $\overline{D}_0(\rho _0)$ the diffusion coefficient renormalised by the interactions. Here, we will assume that the density fluctuations are short enough in time and space, such that $\bar {\rho }(x,y,t)$ may be treated as a constant, and such that the same renormalisation applies here quickly enough to yield a local diffusion constant depending on the local average density $\overline{D}_0(\bar {\rho })$. This assumption is basically the same as that used in the macroscopic fluctuation theory by Bertini et al. (2015), where it is argued that a coarse-grained hydrodynamic description of a system simply leads to an equation of the type (4.1), but with a local diffusion constant and mobility (that should be coupled to the external forces) that depends on the local density field. Here, we assume that we can write a coarse-grained equation for the evolution of the average density $\bar {\rho }$ as

(4.3)\begin{equation} \frac{\partial{\bar{\rho}}}{\partial t} =- \boldsymbol{ \nabla} \boldsymbol{\cdot}\left( - \overline{D}_0(\bar{\rho})\,\boldsymbol{ \nabla} \bar{\rho} - \frac{\overline{D}_0(\bar{\rho})}{k_B T}\,\bar{\rho}\,\boldsymbol{ \nabla} ( \mathcal{V}_{int} * \bar{\rho} ) \right). \end{equation}

As the field $\bar {\rho }(x,y,t)$ is smooth and the interaction is short-range compared to all the other length scales in the system, we can make the local approximation

(4.4)\begin{equation} \mathcal{V}_{int}(x,y) = \frac{E_0}{\rho_0}\,\delta(x)\,\delta(y), \end{equation}

such that the Fokker–Planck equation for the particle density simplifies to

(4.5)\begin{equation} \frac{\partial{\bar{\rho}}}{\partial t} =- \boldsymbol{ \nabla}\boldsymbol{\cdot} \left( - \overline{D}_0(\bar{\rho})\,\boldsymbol{ \nabla} \bar{\rho} - \overline{D}_0(\bar{\rho})\,\frac{E_0}{k_B T}\,\frac{\bar{\rho}}{\rho_0}\,\boldsymbol{ \nabla} \bar{\rho} \right), \end{equation}

which is nonlinear in $\bar {\rho }$.

Our goal is now to obtain an expression for the average density $\bar {\rho }(x,y,t)$. We seek a solution beyond the trivial mean field (where $\bar {\rho }(x,y,t) = \rho _0$), which makes our approach similar in essence to other derivations of particles on lattices (Illien et al. 2018; Rizkallah et al. 2022). Notice that for the height fluctuation profile considered here, the density field is stationary in the frame of reference where the wall is static; i.e. making the change of variables $x' = x - v_{wall}\,t$ and dropping the prime signs yields

(4.6)\begin{equation} 0 = v_{wall}\,\frac{\partial{\bar{\rho}(x,y)}}{\partial x} - \boldsymbol{ \nabla} \boldsymbol{\cdot}\left( - \overline{D}_0(\bar{\rho})\,\boldsymbol{ \nabla} \bar{\rho}(x,y) - \overline{D}_0(\bar{\rho})\,\frac{E_0}{k_B T}\,\frac{\bar{\rho}}{\rho_0}\,\boldsymbol{ \nabla} \bar{\rho} \right). \end{equation}

Integrating vertically, and using exactly the same arguments as in § 2.2.2 for the boundary conditions, shows that

(4.7)\begin{align} J = \int_0^{h(x,t)} \text{d} y \left( - v_{wall}\,\bar{\rho}(x,y) - \overline{D}_0(\bar{\rho})\,\frac{\partial{\bar{\rho}(x,y)}}{\partial x} - \overline{D}_0(\bar{\rho})\,\frac{E_0}{k_B T}\,\frac{\bar{\rho} (x,y)}{\rho_0}\,\frac{\partial{\bar{\rho}(x,y)}}{\partial x} \right) , \end{align}

where $J$, the vertically integrated longitudinal flux, is a constant to be determined.

To solve (4.7), we first assume that we can make the Fick–Jacobs approximation $\bar {\rho }(x,y) \simeq \bar {\rho }(x)$, where the density profile depends only on the longitudinal coordinate, which is correct to first order in $H/L$ and can be demonstrated using lubrication theory. Second, we assume that the average density profile has only small variations around its mean value $\rho _0$, and that the variations are of order $\varepsilon = h_0/H$. We therefore seek a perturbative solution as $\bar {\rho }(x,y) \simeq \rho _0 + \varepsilon \,\rho _1(x) + \cdots$, where $\rho _0 = N/(L H)$ is the mean particle density in the simulation. We can then expand (4.7) in powers of $\varepsilon$. At the lowest order, we find the value of the vertically integrated drift $J = - v_{wall}\,\rho _0 H$. At the next order, we obtain a closed set of equations for the perturbation $\rho _1$,

(4.8)\begin{equation} \left.\begin{gathered} - H\,v_{wall}\,\rho_0 \varepsilon \cos(k_0 x) = H\left[- \overline{D}_0(\rho_0)\left( 1 + \frac{E_0}{k_B T}\right)\frac{\partial{\rho_1(x)}}{\partial x} - v_{wall}\,\rho_1(x) \right] , \\ \int_0^L \rho_1 (x)\,H \,\text{d}x = 0, \\ \rho_1 (x + L) = \rho_1 (x), \end{gathered}\right\} \end{equation}

with relevant periodic boundary conditions, and vanishing integral of $\rho _1$ since the average density should be conserved. We abbreviate $\overline{D}_0 = \overline{D}_0(\rho _0)$ in the following.

Interestingly, we find in (4.8) that the density profile relaxes with an effective diffusion coefficient

(4.9)\begin{equation} D^\star \equiv \overline{D}_0 \left( 1 + \frac{E_0}{k_B T}\right)= \overline{D}_0\,\frac{\chi_T^{id}}{\chi_T^{{\star}}} \end{equation}

that is the collective diffusion coefficient, as it characterises how the fluid of interacting particles relaxes, not how a single particle diffuses. Equation (4.9) is also expected for hard spheres (Hess & Klein 1983; Lahtinen et al. 2001). If interactions are repulsive (i.e. $E_0>0$), then the fluid becomes more incompressible with $\chi _T^\star <\chi _T^{id}$, and density inhomogeneities relax faster than they would in the ideal gas, since $D^\star >\overline{D}_0$. Conversely, attractive interactions (i.e. $E_0<0$) increase compressibility and stabilise density inhomogeneities. This property is absolutely essential to understand the long-time transport properties in a fluctuating medium.

Solving (4.8), we find

(4.10)\begin{equation} \bar{\rho}(x) = \rho_0 + \varepsilon\,\rho_1(x) = \rho_0\left( 1 - \frac{h_0}{H}\,\frac{Pe^{{\star}}}{1 + (Pe^{{\star}})^2}\left[ \sin(k_0x) + Pe^{{\star}} \cos(k_0x) \right] \right), \end{equation}

and a new and natural Péclet number characterising the fluctuations emerges,

(4.11)\begin{equation} Pe^{{\star}} = \frac{v_{wall}}{D^\star k_0} = \frac{\omega_0}{D^\star k_0^2}, \end{equation}

and takes into account the enhanced collective diffusion $D^\star$ due to repulsive interactions. Notice that the solution can be rewritten further as

(4.12)\begin{equation} \bar{\rho}(x) = \rho_0 + \varepsilon\,\rho_1(x) = \rho_0\left( 1 + \frac{h_0}{H}\,\frac{Pe^{{\star}}}{\sqrt{1 + (Pe^{{\star}})^2}} \cos(k_0x + \varphi) \right), \end{equation}

where $\varphi = {\rm \pi}- \arctan ( 1/Pe^{\star } )$. This shows that the perturbation in the density of particles due to the fluctuating interface propagates with a phase $\varphi$ that is characterised by how fast particles relax as a group.

We plot the resulting longitudinal probability density profiles $p(x) = [\rho _0 + \varepsilon \, \rho _1(x)]/L\rho _0$, where $\rho _1(x)$ is given by (4.10), in figures 4(c,d), and find excellent agreement with the numerical results. More specifically, particles accumulate at the constriction at high forcings $v_{wall}$ (figure 4c, light purple) or at low particle densities (figure 4d, yellow and orange). If collective effects are weak, then the time scale for density fluctuations to diffuse across a bulge $1/D^\star k_0^2$ is larger than the time that it takes a bulge to move, $1/v_{wall}\,k_0 = 1/\omega _0$, or equivalently, $Pe^{\star } \gg 1$. Clearly, from (4.12), the phase is $\varphi \simeq {\rm \pi}$ and the interface squeezes particles exactly out of phase, i.e. in the constriction. This can also be understood in the conservation of mass (4.8), where at large $v_{wall}$, we simply have $(H - h(x))\,v_{wall}\,\rho _0 = \varepsilon \,\rho _1(x)\,v_{wall}$ such that $\rho _1(x) \sim 1/h(x)$. Similarly as for the isolated particles in § 2, this is a ‘traffic jam’ effect: in the referential where the wall is fixed, particles trying to move with the fast velocity $-v_{wall}$ accumulate where the road is narrow, i.e. at the constriction.

However, there is a trade-off between accumulation due to wall speed, and increased local repulsion of particles in the accumulation region. We observe that at higher particle densities $\rho _0$, the accumulation effect is tempered (figure 4d, orange to purple). When particle repulsion increases (due to either increased particle density $\rho _0$ or interaction strength $\alpha >0$), the local repulsion dominates, resisting compression of the gas in the narrow constriction. This is coherent with the fact that the compressibility $\chi _T^{\star }/\chi _T^{id}$ decreases with particle density $\rho _0$ (or with the strength of the repulsive interactions $\alpha$). As a result, collective effects are strong, $Pe^{\star } \ll 1$, density fluctuations relax faster than the interface's motion, and the density profile flattens out.

We now turn to understand how this local distribution of particles affects the long-time transport properties of a tracer.

4.2. Effective diffusion and drift with interacting particles

4.2.1. Equation for the diffusion of the tracer

We now consider the motion of a tracer in this gas of soft-core interacting Brownian particles, and we infer how the confined fluctuating environment modifies the long-time tracer motion. Based on reasoning similar to that given above, the probability distribution $\rho _{tr}(x,y,t)$ of finding the tracer particle at position $(x,y)$ at time $t$ satisfies

(4.13) \begin{equation} \frac{\partial{\rho_{tr}(x,y,t)}}{\partial t} =- \boldsymbol{ \nabla} \boldsymbol{\cdot} \left(- \overline{D}_0(\bar{\rho})\,\boldsymbol{\nabla} \rho_{tr} - \overline{D}_0(\bar{\rho})\, \frac{E_0}{k_B T}\,\frac{\rho_{tr}}{\rho_0}\,\boldsymbol{ \nabla} \bar{\rho} \right). \end{equation}

Again looking at the longitudinal transport within the Fick–Jacobs framework, the equation on the marginal distribution $p_{tr}(x,t) = \int _0^{h(x,t)} \rho _{tr}(x,y,t) \,\textrm {d} y$ is given by (2.10):

(4.14) \begin{equation} \frac{\partial{p_{tr}(x,t)}}{\partial t} =- \frac{\partial{}}{\partial x} \left(- \overline{D}_0(\bar{\rho})\,\frac{\partial{p_{tr}}}{\partial x} + \overline{D}_0(\bar{\rho})\,p_{tr}\,\frac{\partial{\ln h(x,t)}}{\partial x} - \overline{D}_0(\bar{\rho})\,\frac{E_0}{k_B T}\,\frac{p_{tr}}{\rho_0}\,\frac{\partial{\bar{\rho}}}{\partial x} \right). \end{equation}

To understand how the tracer's motion is modified by interactions with other particles and by the fluctuating boundary, we expand (4.14) in powers of $\varepsilon =h_0/H$, and obtain the particle's long-time diffusion coefficient and drift, at order $\varepsilon ^2$, as was done in § 2.3. Since the average particle density depends on space, the local diffusion coefficient of the tracer also depends on space as $\overline{D}_0(\bar {\rho }(x,t)) = \overline{D}_0(\rho _0) + \varepsilon \,\overline{D}_0'(\rho _0)\,({\rho _1(x,t)}/{\rho _0}) + O(\varepsilon ^2)$. Here, since the prescribed wall fluctuations $h(x,t)$ are periodic in space, it is possible to obtain implicit expressions for the long-time effective diffusion coefficient $D_{eff}$ and drift $V_{eff}$ using the approach reported in Reimann et al. (2001) and in Guérin & Dean (2015), and resolve their explicit expressions after cumbersome expansions in the small parameter $\varepsilon$ (see § 3 of the supplementary material). In the cases explored here, the local change in diffusion coefficient is small, $\overline{D}_0'(\rho _0)/\overline{D}_0(\rho _0)\,\rho _0 \ll 1$, and we find that its impact on the tracer dynamics can be neglected. We can therefore use the explicit perturbation theory results in § 2.3, assuming $\overline{D}_0(\bar {\rho }(x,t)) \simeq \overline{D}_0(\rho _0) \equiv \overline{D}_0$. As expected, the perturbation theory and the periodic framework approach of Reimann et al. (2001) provide exactly the same results, to leading order in $\varepsilon$. For the sake of completeness, we will nonetheless derive, in a future work, a general perturbation theory with explicit results for Fokker–Planck equations with diffusion and drift with arbitrary dependence on space and time in any dimension.

4.2.2. Local drift of a tracer in a bath of interacting particles

The equation of motion (4.14) can be simplified to

(4.15)\begin{equation} \frac{\partial{p_{tr}(x,t)}}{\partial t} =- \frac{\partial{}}{\partial x} \left(- \overline{D}_0\,\frac{\partial{p_{tr}(x,t)}}{\partial x} + v(x,t)\,p_{tr}(x,t) \right), \end{equation}

where the local longitudinal drift is

(4.16)\begin{equation} v(x,t) = \overline{D}_0\,\frac{1}{H}\,\frac{\partial{h(x,t)}}{\partial x} - \overline{D}_0\,\frac{E_0}{k_B T}\,\frac{h_0}{H}\,\frac{1}{\rho_0}\, \frac{\partial{\rho_1(x,t)}}{\partial x} \end{equation}

at lowest order in $\varepsilon = h_0/H$. We can verify that terms of $O(\varepsilon ^2)$ in the local drift vanish when averaged over one period and hence do not contribute to the renormalization of the long-time transport properties. The last term of (4.16) quantifies the effect of interactions on local tracer motion, and was not present in the ideal gas case. It indicates that particles drift away from accumulation regions because of repulsive interactions. Finally, the magnitude of the effect is proportional to the interaction strength $E_0/k_B T = (\chi _T^{id} - \chi _T^{\star })/\chi _T^{\star }$, or to the relative compressibility of the fluid.

We compute in simulations the mean local particle velocity in the interacting gas of particles; see figure 5. With increasing incompressibility (increasing $\rho _0$ in figures 5a–c), a local drift emerges that carries particles away from the accumulation region located at the channel constriction. As expected with increasing interactions, the velocity field approaches that of the incompressible fluid (reported for comparison in figure 5d, already shown in figure 2d).

Figure 5. Local drift and particle density profiles within the channel for varying mean particle densities. (a–c) Simulation results for increasing mean density for simulations with interacting particles. (d) Simulation results for the incompressible fluid case. Note that (d) corresponds to figure 2(d) and is shown as a reminder. The colour scale for the density profiles is shared between all plots. Yellow (purple) regions indicate regions of high (low) density. The blue arrows correspond to the local velocity of particles (averaged over time), and the length of the arrows is proportional to the magnitude of the velocity. The scale of the arrows is arbitrary but is the same across all plots, hence small arrows in (a) indeed indicate very weak velocity fields compared to (c) or (d). The interaction strength for (a–c) is $\alpha = 1 k_B T/\ell _0^2$, and other numerical parameters are similar to those in figure 3. For all plots, $v_{wall} = 0.5 \ell_0 /\tau_0$.

These results are confirmed by the analytical prediction, using the expression for $\rho _1$ in (4.10). Recalling that the collective Péclet number is $Pe^{\star } = \omega _0^2 / D^\star k_0$, and denoting $\bar {Pe}=\omega _0^2 / \overline{D}_0 k_0$ as the Péclet number of a tracer in medium of density $\rho _0$, we obtain in Fourier space

(4.17)\begin{equation} \tilde{v}(k,\omega) = \frac{\tilde{h}}{H}\,\frac{{\rm i} \overline{D}_0 k (1 + Pe^{{\star}} \bar{Pe}) + (Pe^{{\star}}/\bar{Pe} - 1)\omega/k }{{1 + (Pe^{{\star}})^2}}. \end{equation}

Equation (4.17) interpolates perfectly between the ideal gas and the incompressible fluid. When interactions vanish, $D^\star =\overline{D}_0=D_0$ and $Pe^{\star }=\bar {Pe}=Pe$, which yields the ideal gas expression for the local velocity (2.16), where $\tilde {v}(k,\omega ) = \textrm {i} D_0 k \tilde {h}/H$. Reciprocally, the incompressible fluid limit is $D^\star \to \infty$, which implies $Pe^{\star } \rightarrow 0$, hence ${\tilde {v}(k,\omega ) = (\textrm {i} \overline{D}_0 k - \omega /k)\tilde {h}/H}$, as found already in (2.22) but with $\overline{D}_0$ instead of $D_0$ since the interactions have renormalised the bare diffusion coefficient of the tracer.

4.2.3. Effective long-time diffusion and mean drift

We use the perturbation results in (2.14) and (2.15) to obtain the long-time diffusion coefficient and drift. From (4.17), the spectrum of the local velocity fluctuations is

(4.18)\begin{equation} S(k,\omega) = \frac{\overline{D}_0^2 k^2}{(2{\rm \pi})^2}\,\frac{1 + \bar{Pe}^2}{1 + (Pe^{{\star}})^2}\,\frac{h_0^2}{4 H^2} \left( \delta(k+k_0)\,\delta(\omega - \omega_0) + \delta(k-k_0)\,\delta(\omega + \omega_0) \right), \end{equation}

and the long-time transport properties are

(4.19)\begin{gather} \frac{D_{eff}^{\textit{int. gas}}}{ \overline{D}_0} = 1 - \frac{h_0^2}{2H^2}\,\frac{1 - 3 \bar{Pe}^2}{1 + \bar{Pe}^2}\,\frac{1}{1 + (Pe^{{\star}})^2}, \end{gather}

(4.20)\begin{gather}\frac{V_{eff}^{\textit{int. gas}}}{v_{wall}} = \frac{h_0^2}{2H^2}\,\frac{1}{ 1 + (Pe^{{\star}})^2}. \end{gather}

Equations (4.19) and (4.20) form the main analytic result of the paper and give the transport properties of tracers in a fluid with arbitrary compressibility. We recall that the modified Péclet number is connected to the compressibility $Pe^{\star } = \bar {Pe} \chi _T^{\star }/\chi _T^{id}$, where $\bar {Pe}$ is the Péclet number characterising the fluctuations. For weak interactions, where $\rho _0, \alpha \rightarrow 0$, the compressibility is that of the ideal gas $\chi _T^{\star } \rightarrow \chi _T^{id}$, and $Pe^{\star } \rightarrow Pe$. In this limit, we recover the ideal gas results (2.18) and (2.19). Reciprocally, for strong interactions, $\rho _0, \alpha \rightarrow \infty$, $\chi _T^{\star } \rightarrow 0$, hence $Pe^{\star } \rightarrow 0$ and we recover the incompressible fluid results of (2.24) and (2.25). Note that in the limit $\rho _0\to \infty$, $\overline{D}_0(\rho _0)$ is still finite; see § A.4. Equations (4.19) and (4.20) therefore interpolate between the ideal gas and the incompressible fluid cases.

Despite the mean-field-like assumptions made, our analytic theory summarised in (4.19) and (4.20) perfectly captures the simulation results (see the solid lines in figures 3 and 10), and confirms transport mechanisms in this complex environment. With increasing particle density, particle collisions push particles away from the accumulation region, further into the bulges. Particle–wall collisions then push and disperse particles. In this complex environment, and in contrast with the paradigm where crowded environments slow down diffusion (Lekkerkerker & Dhont 1984; Lowen & Szamel 1993; Dean & Lefèvre 2004), here particle collisions or interactions favour mixing.

5. Discussion and conclusion

In this work, we have explored the impact of crowding on the long-time transport properties of particles in fluctuating channels. Our simulation results show a broad range of behaviours that are well captured by our explicit analytic theory based on a perturbation expansion in the wall fluctuation amplitude $h_0/H$. The results are best described in terms of a Péclet number characterising the fluctuations, $Pe= \omega _0/D_0 k_0^2 = 2{\rm \pi} v_{wall}/D_0 L$. At low $Pe \ll 1$, corresponding to fixed interfaces, all fluids behave similarly: particles are slowed down by constrictions, and the effective diffusion is decreased. At intermediate $Pe \simeq 1$, particle–wall collisions stir particles, and the effective diffusion is increased. This effect persists until the wall moves so fast that particles no longer have time to reach the moving bulges and accumulate in the constrictions. In this final regime, $Pe \gg 1$, the effective diffusion is unchanged. The accumulation regime arises for higher and higher $Pe$ values for increasing particle–particle interactions, i.e. for increasing incompressibility that resists accumulation.

5.1. Collisions enhance diffusion

One of the main findings of our work is that here, both numerically and analytically, we have demonstrated that increasing repulsive interactions or collisions between particles can enhance the late-time diffusion coefficient and the mean drift characterising the dispersion of a tracer particle in fluctuating channels. This is in contrast with the intuition that collisions in equilibrium reduce the diffusion coefficient (Lekkerkerker & Dhont 1984; Lowen & Szamel 1993). The mechanism of diffusion enhancement is in fact rather simple: collisions or repulsive interaction help to push particles closer to the walls. Eventually, wall–particle collisions help mixing. Since this mechanism is rather straightforward, we expect it to be broadly applicable – for example, beyond lubrication approximation or for hard-core repulsive interactions. Such effects could be explored within our framework or alternatively using dynamic density functional theory (Marconi & Tarazona 1999). In more detailed physical settings, such as with hard-core interactions, other effects would likely also come into play; for example, the accessible volume in the channel is smaller for larger particles (Riefler et al. 2010; Suárez et al. 2015), and hydrodynamic effects become important (Yang et al. 2017).

We now put our results in a broader context. In the Introduction, we recalled that diffusion of a driven, out-of-equilibrium tracer in a bath of interacting particles is enhanced by repulsive interactions (Bénichou et al. 2013, 2018; Démery et al. 2014; Illien et al. 2018). In a confined channel, thermal fluctuations of the wall could possibly enhance the diffusion coefficient of particles, as we have seen in the limiting case of an incompressible fluid (Marbach et al. 2018). Another recent work finds that diffusion of odd-diffusing interacting particles is enhanced with increasing densities (Kalz et al. 2022). While the physical set-up in Kalz et al. (2022) is very different from ours, the mathematical similarities that lead to diffusion enhancement are striking (for example, comparing their Eq. (9) with our (2.24)), and one might speculate that there exists a universal framework to understand these effects under the same light.

5.2. Open questions on fluctuating channels

Beyond the question asked here, namely of understanding how crowding affects transport in fluctuating channels, there are many open fundamental questions. For example, boundaries are not necessarily repulsive and smooth. Surface rugosity would lead to Taylor dispersion in an incompressible fluid (Marbach & Alim 2019; Kalinay 2020), but how would surface rugosity of the wall potential induce Taylor dispersion in the fluid of interacting particles? Attraction at the boundaries also leads to surprising speed-up of diffusion in corrugated static channels (Alexandre et al. 2022). Is this speed-up enhanced further by fluctuations? Considering more complex fluids would likely modify the distribution of particles along the channel, and hence the transport properties. For example, normal stresses close to the boundaries in odd viscous fluids bring particles close to the walls (Hargus et al. 2020; Fruchart, Scheibner & Vitelli 2023). Could this phenomenon enhance dispersion? Down the scales, molecular (Yoshida et al. 2018) or quantum (Kavokine, Bocquet & Bocquet 2022; Coquinot, Bocquet & Kavokine 2023) effects enhance the mobility of individual molecules; how would these effects combine with mechanical fluctuations? With the advent of highly sensitive techniques to probe the motion of particle near surfaces in soft and increasingly complex environments (Zhang et al. 2020; Sarfati et al. 2021; Vilquin et al. 2022), one might expect to answer some of these questions in the light of further experimental results.