## 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 Reference Bénichou, Illien, Oshanin, Sarracino and Voituriez2019; Tomkins, Hughes & Morris Reference Tomkins, Hughes and Morris2021), to improved filtration of liquids or gases, including modern desalination techniques and oil recovery (Werber, Osuji & Elimelech Reference Werber, Osuji and Elimelech2016; Marbach & Bocquet Reference Marbach and Bocquet2019). 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 Reference Kavokine, Netz and Bocquet2021). 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 Reference Zwanzig1992). This effect is well captured by the so-called Fick–Jacobs formalism (Jacobs Reference Jacobs1935; Reguera & Rubí Reference Reguera and Rubí2001; Kalinay & Percus Reference Kalinay and Percus2006; Rubi Reference Rubi2019), 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 Reference Dagdug, García-Chung and Chacón-Acosta2016; Chávez, Chacón-Acosta & Dagdug Reference Chávez, Chacón-Acosta and Dagdug2018). Recent work suggests that the approach is also adapted to study fluid flow through biological membranes (Arango-Restrepo & Rubi Reference Arango-Restrepo and Rubi2020). Such entropic contributions induce significant corrections to transport in microfluidic (Yang *et al.* Reference Yang, Liu, Li, Marchesoni, Hänggi and Zhang2017) or biological (Rubí *et al.* Reference Rubí, Lervik, Bedeaux and Kjelstrup2017) 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 Reference Leroy, Rousseau and Fuchs2004; Haldoupis *et al.* Reference Haldoupis, Watanabe, Nair and Sholl2012), and more recently, enhanced water diffusion in solid state pores via phonon-fluid coupling (Ma *et al.* Reference Ma, Grey, Shen, Urbakh, Wu, Liu, Liu and Zheng2015, Reference Ma, Tocci, Michaelides and Aeppli2016; Cao, Wang & Ma Reference Cao, Wang and Ma2019; Noh & Aluru Reference Noh and Aluru2021). 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 Reference Marbach, Dean and Bocquet2018). 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 Reference Taylor1953; Aris Reference Aris1956; Hydon & Pedley Reference Hydon and Pedley1993). 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 Reference Chakrabarti and Saintillan2020) and in experiments in fluctuating porous matrices (Sarfati, Calderon & Schwartz Reference Sarfati, Calderon and Schwartz2021). Such enhancement of dispersion properties is relevant in a number of biological contexts, such as in blood vessels (Masri, Puelz & Riviere Reference Masri, Puelz and Riviere2021), slime mould vasculature (Marbach & Alim Reference Marbach and Alim2019), the gut (Codutti, Cremer & Alim Reference Codutti, Cremer and Alim2022), and near molecular motors (Evans, Krause & Feringa Reference Evans, Krause and Feringa2021), and are of general relevance in plants (Tomkins *et al.* Reference Tomkins, Hughes and Morris2021).

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 Reference Lekkerkerker and Dhont1984; Lowen & Szamel Reference Lowen and Szamel1993; Dean & Lefèvre Reference Dean and Lefèvre2004). However, subtle effects may arise if the crowded medium is driven out of equilibrium (Zaccone *et al.* Reference Zaccone, Dorsaz, Piazza, De Michele, Morbidelli and Foffi2011; Lappala, Zaccone & Terentjev Reference Lappala, Zaccone and Terentjev2013). 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.* Reference Bénichou, Illien, Oshanin and Voituriez2013, Reference Bertini, De Sole, Gabrielli, Jona-Lasinio and Landim2018; Démery, Bénichou & Jacquin Reference Démery, Bénichou and Jacquin2014; Illien *et al.* Reference Illien, Bénichou, Oshanin, Sarracino and Voituriez2018). 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.* Reference Illien, Bénichou, Oshanin, Sarracino and Voituriez2018). Such effects may be exploited for active microrheology within spatially corrugated channels (Puertas, Malgaretti & Pagonabarraga Reference Puertas, Malgaretti and Pagonabarraga2018; Malgaretti, Puertas & Pagonabarraga Reference Malgaretti, Puertas and Pagonabarraga2022).

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.* Reference Obliger, Valdenaire, Ulm, Pellenq and Leyssale2019), sometimes even exhibiting a maximum (Pireddu *et al.* Reference Pireddu, Pazzona, Demontis and Załuska-Kotur2019) 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.* Reference Obliger, Bousige, Coasne and Leyssale2023). 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 1*a*), 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.

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.* Reference Marbach, Dean and Bocquet2018). 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 1*b*). 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.* Reference Marbach, Dean and Bocquet2018), 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.* Reference Sznitman, Heimsch, Heimsch, Rusch and Rösgen2007; Dong, Yang & Zhu Reference Dong, Yang and Zhu2022). 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.* Reference Marbach, Dean and Bocquet2018).

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,

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.* Reference Marbach, Dean and Bocquet2018). 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$.

Interestingly, the long-time diffusion coefficient $D_{eff}$ exhibits a non-monotonic dependence on the Péclet number $Pe$ (yellow dots in figure 2*a*), 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 Reference Jacobs1935; Zwanzig Reference Zwanzig1992; Reguera & Rubí Reference Reguera and Rubí2001; Kalinay & Percus Reference Kalinay and Percus2006; Burada *et al.* Reference Burada, Schmid, Reguera, Rubi and Hänggi2007; Mangeat, Guérin & Dean Reference Mangeat, Guérin and Dean2017; Marbach *et al.* Reference Marbach, Dean and Bocquet2018; Rubi Reference Rubi2019). 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 2*c*). 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.* (Reference Marbach, Dean and Bocquet2018), 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

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

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 Reference Chaikin, Lubensky and Witten1995) or peristaltic motion of vasculature (Alim *et al.* Reference Alim, Amselem, Peaudecerf, Brenner and Pringle2013; Marbach & Alim Reference Marbach and Alim2019). 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

where we used the simple and useful relation

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

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

explicitly, and using (2.6) again, we find

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)$:

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 (Reference Mercer and Roberts1990, Reference Mercer and Roberts1994) and Marbach & Alim (Reference Marbach and Alim2019). 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 Reference Jacobs1935; Zwanzig Reference Zwanzig1992; Reguera & Rubí Reference Reguera and 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.* Reference Marbach, Dean and Bocquet2018).

#### 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

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

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

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)

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.* (Reference Marbach, Dean and Bocquet2018), 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$,

This yields a spectrum for the local drift as

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

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 Reference Zwanzig1992; Jacobs Reference Jacobs1935; Reguera & Rubí Reference Reguera and Rubí2001; Marbach *et al.* Reference Marbach, Dean and Bocquet2018; Rubi Reference Rubi2019). 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.* Reference Marbach, Dean and Bocquet2018). 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

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 Reference Marbach and Alim2019; Chakrabarti & Saintillan Reference Chakrabarti and Saintillan2020). 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

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

We then obtain

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.* (Reference Marbach, Dean and Bocquet2018).

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 2*d*), 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 1*b*), 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,

where $r$ is the inter-particle distance, $d_c$ is a critical distance characterising the radius of the interaction (see figure 1*b*), 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.* Reference Bénichou, Illien, Oshanin and Voituriez2013; Suárez, Hoyuelos & Mártin Reference Suárez, Hoyuelos and Mártin2015) 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 Reference Pàmies, Cacciuto and Frenkel2009; Démery *et al.* Reference Démery, Bénichou and Jacquin2014; Antonov, Ryabov & Maass Reference Antonov, Ryabov and Maass2021). 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 Reference Hansen and McDonald2013):

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

We calculate

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

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.* Reference Démery, Bénichou and Jacquin2014) (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.* Reference Démery, Bénichou and Jacquin2014; 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.

### 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 4*a*,*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 2*d*). However, at higher wall velocities (figure 4*b*), 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.* Reference Suárez, Hoyuelos and Mártin2015; Suárez, Hoyuelos & Mártin Reference Suárez, Hoyuelos and Mártin2016), which share a similar geometry in the reference frame where the wall is static.

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 4*c*). It flattens out with increasing particle density $\rho _0$ (see figure 4*d*), 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 Reference Dean1996; Kawasaki Reference Kawasaki1998)

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:

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.* (Reference Démery, Bénichou and Jacquin2014), 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.* (Reference Bhattacharjee and Datta2015), 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

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

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

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.* Reference Illien, Bénichou, Oshanin, Sarracino and Voituriez2018; Rizkallah *et al.* Reference Rizkallah, Sarracino, Bénichou and Illien2022). 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

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

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$,

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

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 Reference Hess and Klein1983; Lahtinen *et al.* Reference Lahtinen, Hjelt, Ala-Nissila and Chvoj2001). 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

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

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

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 4*c*, light purple) or at low particle densities (figure 4*d*, 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 4*d*, 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

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):

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.* (Reference Reimann, Van den Broeck, Linke, Hänggi, Rubi and Pérez-Madrid2001) and in Guérin & Dean (Reference Guérin and Dean2015), 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.* (Reference Reimann, Van den Broeck, Linke, Hänggi, Rubi and Pérez-Madrid2001) 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

where the local longitudinal drift is

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 5*a*–*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 5*d*, already shown in figure 2*d*).

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

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

and the long-time transport properties are

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 Reference Lekkerkerker and Dhont1984; Lowen & Szamel Reference Lowen and Szamel1993; Dean & Lefèvre Reference Dean and Lefèvre2004), 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 Reference Lekkerkerker and Dhont1984; Lowen & Szamel Reference Lowen and Szamel1993). 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 Reference Marconi and Tarazona1999). 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.* Reference Riefler, Schmid, Burada and Hänggi2010; Suárez *et al.* Reference Suárez, Hoyuelos and Mártin2015), and hydrodynamic effects become important (Yang *et al.* Reference Yang, Liu, Li, Marchesoni, Hänggi and Zhang2017).

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.* Reference Bénichou, Illien, Oshanin and Voituriez2013, Reference Bertini, De Sole, Gabrielli, Jona-Lasinio and Landim2018; Démery *et al.* Reference Démery, Bénichou and Jacquin2014; Illien *et al.* Reference Illien, Bénichou, Oshanin, Sarracino and Voituriez2018). 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.* Reference Marbach, Dean and Bocquet2018). Another recent work finds that diffusion of odd-diffusing interacting particles is enhanced with increasing densities (Kalz *et al.* Reference Kalz, Vuijk, Abdoli, Sommer, Löwen and Sharma2022). While the physical set-up in Kalz *et al.* (Reference Kalz, Vuijk, Abdoli, Sommer, Löwen and Sharma2022) 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 Reference Marbach and Alim2019; Kalinay Reference Kalinay2020), 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.* Reference Alexandre, Mangeat, Guérin and Dean2022). 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.* Reference Hargus, Klymko, Epstein and Mandadapu2020; Fruchart, Scheibner & Vitelli Reference Fruchart, Scheibner and Vitelli2023). Could this phenomenon enhance dispersion? Down the scales, molecular (Yoshida *et al.* Reference Yoshida, Kaiser, Rotenberg and Bocquet2018) or quantum (Kavokine, Bocquet & Bocquet Reference Kavokine, Bocquet and Bocquet2022; Coquinot, Bocquet & Kavokine Reference Coquinot, Bocquet and Kavokine2023) 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.* Reference Zhang, Bertin, Arshad, Raphael, Salez and Maali2020; Sarfati *et al.* Reference Sarfati, Calderon and Schwartz2021; Vilquin *et al.* Reference Vilquin, Bertin, Raphaël, Dean, Salez and Mcgraw2022), one might expect to answer some of these questions in the light of further experimental results.

## Supplementary material

Supplementary material is available at https://doi.org/10.1017/jfm.2023.640

## Acknowledgements

The authors are grateful for fruitful discussions with L. Bocquet, J.-P. Hansen, P. Levitz, M.-T. Hoang Ngoc, G. Pireddu, B. Rotenberg, B. Sprinkle and A. Thorneywork. S.M., R.Z. and Y.W. thank the Applied Math Summer Undergraduate Research Experience Program of the Courant Institute for putting them in contact. S.M. would like to thank the Institut d’Études Scientifiques de Cargése for hosting the Transport in Narrow Channels workshop that led to inspiring discussions for this work. R.Z. would like to thank the Laboratoire MSC Paris and the Center for Data Science ENS Paris for hospitality.

## Funding

This work was supported in part by the MRSEC Program of the National Science Foundation under award number DMR-1420073. S.M. was also supported by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement 839225 MolecularControl. R.Z. was also supported by grant no. NSF DMR-1710163.

## Declaration of interests

The authors report no conflict of interest.

## Data availability statement

The data that support the findings of this study are available upon reasonable request to the author.

## Appendix A. Simulation details

#### A.1. Simulation algorithm

All simulations are performed using a forward Euler stochastic scheme to discretise the overdamped Langevin dynamics of the particles. For a particle $i$ at position ${\boldsymbol {X}_i(t) = (x_i(t),y_i(t))}$, the following position at time $t + \Delta t$ is computed as

where $\boldsymbol {U}_i = \boldsymbol {U} (x_i(t),y_i(t))$ is the background flow field (which is non-zero only in the incompressible case), $\boldsymbol {F}_i = \boldsymbol {F}_{wall} + \sum _{j\neq i} \boldsymbol {F}_{{int}, ij}$ is the sum of the forces exerted by the channel walls, $\boldsymbol {F}_{wall}$, on the particle and by the neighbouring particles $\boldsymbol {F}_{{int}, ij}$ when interactions are present, and $\boldsymbol {G}_i(t)$ is a vector of two independent random numbers drawn from a Gaussian distribution of mean 0 and variance 1. Note that our simulations have been tested with a time step twice as small and yielded no significant difference.

Particles are confined in the channel by means of a potential that exerts a force on the particles only if they reach the boundaries. More precisely, the force exerted by the wall on a particle with coordinate $(x,y)$ is given by $\boldsymbol {F}_{wall}=-\boldsymbol { \nabla } \mathcal {V}_{wall}(x,y)$, with

with $\lambda >0$ a stiffness coefficient (with dimensions of energy over length to the power 4), and where the boundary equations are given by

where $H$ and $h_0$ represent the average height and the variation amplitude of the channel height, respectively. Using a soft confining potential to model the wall is convenient for simulation purposes, as it avoids dealing with reflecting Brownian walks, which carries some challenges (Scala, Voigtmann & De Michele Reference Scala, Voigtmann and De Michele2007). It also allows one to keep a rather large integration time step. Finally, this choice of boundary conditions bears a physical meaning since it considers that the wall is likely made of a slightly penetrable material that would interact with the particles/colloids, and is, therefore, more general than a hard-reflecting boundary condition. Such boundary models have been used extensively in theoretical active matter systems (Solon *et al.* Reference Solon, Fily, Baskaran, Cates, Kafri, Kardar and Tailleur2015; Zakine *et al.* Reference Zakine, Zhao, Knežević, Daerr, Kafri, Tailleur and van Wijland2020; Ben Dor *et al.* Reference Ben Dor, Ro, Kafri, Kardar and Tailleur2022), and as our theory relies only on the presence of a boundary, the results are largely unaffected by a change of potential, as long as the boundary layer of particles at the wall and subjected to the potential repulsion is small compared to all other relevant length scales in the system. The fact that our numerical results (and especially the two-dimensional density profiles) are consistently in agreement with hard-reflecting boundary conditions taken in the analytical computation shows that this choice of boundary conditions is fairly equivalent to hard-reflecting boundaries.

In § 3, we have performed simulations with pairwise interacting particles. The force on particle $i$ exerted by particle $j$ is $\boldsymbol {F}_{{int}, ij} = - \boldsymbol { \nabla } \mathcal {V}_{int}(r_{ij})$, where $r_{ij}$ is the distance between $i$ and $j$. We chose simple repulsive interactions described by the soft potential (3.1), which we recall here as $\mathcal {V}_{int}(r) = \alpha (r - d_c)^2\,\varTheta (r < d_c)$, where $\alpha$ is the interaction strength in units of a spring constant, and $d_c$ is the typical particle diameter (see figure 6). For this quadratic interaction $\mathcal {V}_{int}$, the expression of its vertically integrated version is given by

This expression is used to compute the density profiles $\rho (x)$, but we always use the radial potential $\mathcal {V}_{int}$ in the Monte Carlo simulations.

The system size along the $x$ direction is $L$, and we take periodic boundary conditions in the $x$ direction. We keep track of the full dynamics of the $N$ particles with their absolute position $\boldsymbol {X}_{true}$ to compute their mean-square displacement and mean drift, but interactions are computed using the folded positions in the periodic domain $[0,L)\times (-\infty, \infty )$. The mean density of particles is $\rho _0 = N/(L H)$. To access higher particle densities, we therefore change the total number of particles $N$ in the simulation (see table 1).

We initialise systems from a uniform distribution of particles in the bottom part of the channel ($0< y < H- h_0$). We first let the system equilibrate for a time $t_{eq}\sim H/D_0\sim 100 \tau_0$ within a fixed undulated channel ($v_{wall} = 0$). Then the actual simulation starts with a positive value of $v_{wall}$.

#### A.2. Simulation parameters

The domain characteristics are $L=200\ell _0$, $H=12\ell _0$, $h_0=3\ell _0$ for all the simulations. To simulate non-interacting particles, we take wall stiffness $\lambda =300 k_B T/\ell _0^4$, and the integration time step is $\Delta t=4\times 10^{-3} \tau _0$. For a typical simulation, the total number of iterations is $N_{iter}=5\times 10^6$, and the number of particles is $N=10^5$. To simulate interacting particles, we take $\lambda =10 k_B T/\ell _0^4$, and the integration time step is $\Delta t=4\times 10^{-3} \tau _0$. Typical values of the number of particles $N$ and total number of iterations $N_{iter}$ are shown in table 1.

#### A.3. Simulation analysis

We perform multiple independent simulations to increase statistical resolution. For each value of $v_{wall}$, we perform 10 independent simulations starting from different initial configurations (and different seeds). Symbols in all graphs represent the mean value of the observable, and error bars correspond to one standard deviation over these 10 independent measurements. The mean drift is calculated simply as $V_{eff} = ({1}/{N})\sum _{i=1}^N \langle {(x_{{true}, i} (t) - x_{true}(0))}/{t} \rangle _t$. The mean-square displacement of particles is calculated as ${MSD}(t) = ({1}/{N})\sum _{i=1}^N \langle (x_{{true}, i} (t + t_0) - x_{true}(t_0))^2 \rangle _{t_0}$, where the average is done over initial times $t_0$. The self-diffusion coefficient is then obtained as a least squares linear fit of the mean-square displacement. The parameters that we choose allow us to neglect the finite size corrections due to periodic boundary conditions. Indeed, such corrections scale as ${\sim d_c/L=0.005}$ (Dünweg & Kremer Reference Dünweg and Kremer1991, Reference Dünweg and Kremer1993) or as $\sim (H/L)^2=0.004$ (Simonnin *et al.* Reference Simonnin, Noetinger, Nieto-Draghi, Marry and Rotenberg2017), thus are negligible in the numerical measurements performed here.

#### A.4. Simulation calibration: self-diffusion coefficient $\overline{D}_0(\rho _0)$ of soft-core interacting particles

To calibrate our model, we calculate the long-time self-diffusion coefficient of particles in a fluid of soft-core interacting particles. In figure 7, we report the results of the computation of the diffusion coefficient $\overline{D}_0(\rho _0)$ for fixed flat walls according to different parameters of the interaction: the particle density $\rho _0$, and the interaction strength $\alpha$.

As a self-consistent check, we also compare our simulations to the analytic mean-field theory of Démery *et al.* (Reference Démery, Bénichou and Jacquin2014), and find good agreement between simulation and theory. We recall briefly the analytic formula here. The correction to the bare diffusion coefficient is given by

where $\rho _0=N/V$ is the particle density, $d$ is the number of spatial dimensions, and $\tilde{\mathcal {V}}_{int}(\boldsymbol {k})$ is the Fourier transform of the interaction pair potential, here

Note that at large densities $\rho _0$ or interaction strengths, (A7) is no longer valid, and it has been argued that it should be approached by (Dean & Lefèvre Reference Dean and Lefèvre2004)

In our set-up for $d=2$ and $k=|\boldsymbol {k}|$, we have

where $J_\nu (z)$ is the special Bessel function, and $H_\nu (z)$ is the Struve function

We plot (A7) in figure 7, and compare it to our numerical results. Unsurprisingly, the mean-field approximation starts to fail as the interaction strength $\alpha$ increases and prevents particles from overlapping. To present the results of the long-time diffusion coefficients $D_{eff}$ in confined wiggling spaces relative to $D_0$ (figures 3 and 10), we always use the numerically calculated diffusion coefficient $\overline{D}_0(\rho _0)$ in the flat fixed space.

In the main text, we investigate limit behaviours as $\rho _0 \rightarrow \infty$ and $\alpha \rightarrow \infty$. A convenient Gaussian interaction potential can be considered to gain analytical insights into the diffusion. With $\mathcal {V}_{int}(\boldsymbol {r})=\alpha \ \textrm {e}^{-\boldsymbol {r}^2/(2d_c^2)}$, we have $\tilde {\mathcal {V}}_{int}(\boldsymbol {k})=2{\rm \pi} \alpha d_c^2 \ \textrm {e}^{-d_c^2\boldsymbol {k}^2}$in dimension $d=2$, for instance. The integral in (A9) can thus be approximated when integrating only in the range ${\rho _0\,\tilde {\mathcal {V}}_{int}(\boldsymbol {k})}/{(k_B T)}\gg 1$. The interaction-and-density-dependent cutoff is given by $\varLambda (\rho _0)\sim ({1}/{d_c})\sqrt {\ln (\alpha d_c^2 \rho _0/k_B T)}$, and the diffusion is thus given by

which yields, first, $\overline{D}_0(\rho _0)\to D_0$ for $\rho _0\to \infty$ and $\alpha$ fixed (as expected since the potential landscape becomes flat), and second, $\overline{D}_0(\rho _0)\to 0$ for $\alpha \to \infty$ and $\rho _0$ fixed (as expected for jamming).

#### A.5. Simulation calibration: confinement with the soft wall potential

For each simulation type, we check the penetration depth of the particles in the confining wall. With increasing repulsive interactions (with increased $\rho _0$ or $\alpha$), and since the wall is ‘soft’, particles may be squeezed into the confining soft wall. We then estimate the penetration depth of each fluid within the confining wall by looking at the probability density at the centre of the channel where the constriction is (in the frame of reference where the wall is fixed); see figure 8. We find that indeed, with increasing interactions the penetration depth $\delta h = (0.2- 1) \ell _0$ increases. We record the penetration depth $\delta h$ from figure 8 for each set of numerical parameters, corresponding to the depth for which the probability density is half of its bulk value (dashed black horizontal line). We then use $H = H + 2 \delta h$ (since there is an upper wall and a bottom wall) in all analytical formulas.

## Appendix B. Additional data for the ideal gas

To test the validity of the analytic derivation in § 2 for the transport of isolated particles, we explore here a broader range of simulation parameters. In particular, we go beyond the lubrication approximation and investigate systems for which $H/L \simeq 0.02$ up to $H/L \simeq 0.1$. We report the measured long-time diffusion coefficients $D_{eff}$ and mean drift $V_{eff}$ in figure 9, along with representative plots of the density profile in the frame of reference where the channel wall is fixed. We find that as the width $H$ increases, a $y$ dependence of the stationary density emerges (see figure 9*e*). In fact, diffusion across the channel width can no longer be considered fast with respect to diffusion along length $x$. This corresponds to the progressive breakdown of the lubrication approximation.

Surprisingly, the analytic formulas (2.18) and (2.19) are still in remarkable agreement with simulations, up to $H/L \simeq 0.1$. Slight deviations may be observed for $H/L \simeq 0.1$ (corresponding to $H = 20 \times \ell _0$) on the mean drift, where $V_{eff}/v_{wall} > 0$ even at large $Pe$. This is due to accumulation of particles in the upstream bulge, as they collide and leave a wake of particles, instead of having the time to distribute vertically. As a result, some particles are still in the bulge even at large $Pe$, and are therefore carried, which produces a net mean drift. Interestingly, at very small $H/L \simeq 0.02$ (corresponding to $H = 4 \times \ell _0$), we observe slight deviation from the theory this time at small $Pe$. This is due to the fact that for such systems, the penetration depth in the wall, $\delta h \simeq 0.2 \ell _0$, becomes more and more comparable with the channel height $H$. The effective vertical accessible space $H$ is therefore larger, and the value in (2.18) and (2.19) should be modified appropriately (see § A.5).

## Appendix C. Perturbation theory to obtain long-time transport coefficients

Here, we perform the perturbation theory to obtain the long-time transport coefficients $D_{eff}$ and ${V_{eff}}$ of the general Fokker–Planck equation (2.12).

In the following, it will be easier to work in Fourier space, and we therefore define, for any function $f(x,t)$, the Fourier transform

where the $\int$ sign encompasses integration over space and time. Conversely, the reverse Fourier transform is given by

Performing a Fourier transform on (2.12) yields

which can be written, using the notation $1/(D_0 k^2 + \textrm {i} \omega )\equiv \tilde {p}_0(k,\omega )$, as

We would like to eventually simplify (C4) in a form where diffusion and advection coefficients can be read easily, as suggested by the Fourier transform of (2.11) that yields

The lubrication approximation enables us to simplify further the effective equation on $\tilde {p}(k,\omega )$. Indeed, we consider that the fluctuations at the channel boundary are a perturbation, with small relative amplitude $\varepsilon = h_0/H$. Hence the last term of (C4), containing the advection $v$, can be considered as a perturbation. For example, in the case of the ideal gas,

where we denote $\tilde {h}(k,\omega )$ the Fourier transform of the non-constant part of $h(x,t)$, i.e. $h(x,t)- H$. We therefore seek a solution to (C4) as an expansion in $\varepsilon$, namely, $\tilde {p}(k,\omega ) = \tilde {p}_0(k,\omega ) + \varepsilon \, \tilde {p}_1(k,\omega ) + \varepsilon ^2\,\tilde {p}_2(k,\omega ) +\cdots$.

Additionally, since we seek the behaviour of the solution at long times, it is natural to calculate the noise-averaged solution $\langle \tilde {p}(k,\omega ) \rangle$ where $\langle \cdot \rangle$ denotes the usual average over realisations of the noise. Note that we can also treat the propagating wave case, which is deterministic, in terms of a random field. Consider that we define $h(x,t) = H + h_0 \cos (k_0 x - \omega _0 t + \theta )$, where $\theta$ is a random phase distributed uniformly on $[0,2{\rm \pi} ]$. Clearly, the value of $\theta$ cannot affect the result at late times as it just fixes the initial configuration of the height at time $t = 0$ when the advection diffusion process starts. This choice of $\theta$ is also equivalent, for instance, to choosing an arbitrary initial time $t_0 = \theta /\omega _0$. We thus define $\langle {\cdot } \rangle$ here as a uniform average over $\theta$ on $[0,2{\rm \pi} ]$. Using this convention, we see immediately that $\langle h \rangle = H$. Here, we consider additionally Gaussian fluctuations with mean 0, i.e. the first two moments of the noise completely specify the problem. Thus $\langle \tilde {v}(k,\omega )\rangle = 0$ and

where $S(k,\omega )$ corresponds to the spectrum of the fluctuating advection, as defined in (2.13). Of course, $\langle \tilde {p}_0(k, \omega ) \rangle = \tilde p_0(k, \omega )$.

We now solve iteratively for $\tilde {p}(k,\omega )$. At first order in $\varepsilon$, we obtain

and $\langle \tilde {p}_1(k,\omega ) \rangle = 0$. At second order in $\varepsilon$, we have

and averaging over noise gives

We observe that we can define $\varSigma (k,\omega )$ such that

and $\langle \tilde {p}_2(k,\omega ) \rangle = \tilde {p}_0(k,\omega )\,\varSigma (k,\omega )\, \tilde {p}_0(k,\omega )$.

Pursuing the derivation, one can show that the full solution is a geometric series $\langle \tilde {p}\rangle = \tilde {p}_0 + \varepsilon ^2\,\tilde {p}_0 \varSigma \tilde {p}_0 + \varepsilon ^4\, \tilde {p}_0 \varSigma \tilde {p}_0 \varSigma \tilde {p}_0 +\cdots$ that can be re-summed to obtain

Additional steps may be found in the supplementary information of Marbach *et al.* (Reference Marbach, Dean and Bocquet2018).

#### C.1. Long-time results

From the target long-time expression (C5) and from the re-summed propagator (C12), one can read off the effective long-time diffusion coefficient as

and the mean drift as

Injecting $\varSigma (k,\omega )$ from (C11) into (C13) and (C14) (and dropping the $'$ in the integrals for simplicity), we obtain

and

Now, assuming that the problem is translationally invariant in space and time (which is reasonable considering that the channel is assumed to be infinitely long, and that such an assumption still allows one to model many different situations), this implies that the height correlations satisfy $\langle h(x,t)\,h(x',t') \rangle = C(|x-x'|,|t-t'|)$, hence $S(k,\omega ) = S(-k,-\omega )$. In addition, as the correlation function $C$ is real, the conjugate is $S^*(k,\omega ) = S(-k,-\omega )$. Using time and space invariance, we obtain that also $S$ is real. Using the fact that $S$ is even with respect to both its variables, and plugging back the expression of $\varepsilon = h_0/H$, (C15) simplifies to (2.14), and (C16) to (2.15), of the main text.

## Appendix D. Additional data for a fluid of interacting particles

We present here additional data for interacting particles. We inspect various values of the interaction strength $\alpha$. We present in figure 10 the results of the normalised effective diffusion $D_{eff}/\overline{D}_0(\rho _0)$ and the normalised mean drift ${V_{eff}} /v_{wall}$. We find very similar results whether probing increasing interaction strength $\alpha$ or probing increasing particle density $\rho _0$ (see figure 3). In fact, increasing $\alpha$ also increases the energy scale $E_0(\alpha,\rho _0)$ contained in a volume $d_c^d$ (with $d$ the dimension of space). Hence, since the theory depends mainly on the value of $E_0= ({\rm \pi} /6) \alpha d^4_c \rho_0$, similar effects are observed naturally when increasing $\alpha$ or $\rho_0$.