# How stickiness can speed up diffusion in confined systems

A. Alexandre  Laboratoire Ondes et Matière d’Aquitaine, CNRS/University of
Bordeaux, F-33400 Talence, France  M. Mangeat  Center for Biophysics and
Department for Theoretical Physics, Saarland University, D-66123 Saarbrücken,
Germany  T. Guérin  Laboratoire Ondes et Matière d’Aquitaine, CNRS/University
of Bordeaux, F-33400 Talence, France  D.S. Dean  Laboratoire Ondes et Matière
d’Aquitaine, CNRS/University of Bordeaux, F-33400 Talence, France  Team MONC,
INRIA Bordeaux Sud Ouest, CNRS UMR 5251, Bordeaux INP, Univ. Bordeaux,
F-33400, Talence, France

###### Abstract

The paradigmatic model for heterogeneous media used in diffusion studies is
built from reflecting obstacles and surfaces. It is well known that the
crowding effect produced by these reflecting surfaces slows the dispersion of
Brownian tracers. Here, using a general adsorption desorption model with
surface diffusion, we show analytically that making surfaces or obstacles
attractive can accelerate dispersion. In particular, we show that this
enhancement of diffusion can exist even when the surface diffusion constant is
smaller than that in the bulk. Even more remarkably, this enhancement effect
occurs when the effective diffusion constant, when restricted to surfaces
only, is lower than the effective diffusivity with purely reflecting
boundaries. We give analytical formulas for this intriguing effect in periodic
arrays of spheres as well as undulating micro-channels. Our results are
confirmed by numerical calculations and Monte Carlo simulations.

Determining the transport properties of tracer particles in heterogeneous
media at large time and length scales has applications in a range of physical
problems including fluid mechanics, hydrology, chemical engineering, soft
matter and solid state physics [1, 2, 3, 4]. The effective diffusivity is a
crucial input for problems of mixing [5, 6, 7], sorting [8], chemical delivery
[2, 9] as well as chemical reactions [6, 4]. Spatial variations of diffusion
and advection can lead to drastic changes in the effective diffusion constant
with respect to homogenous systems. Classic examples include Taylor dispersion
in hydrodynamic flows [10], and the decrease of dispersion due to the energy
barriers created by time independent potentials [11]. In a number of systems,
such as porous media [12], zeolites, and biological channels [13], the motion
of tracer particles is hindered by hard (technically speaking reflecting)
boundaries or obstacles. The effective trapping of the tracer in this case is
of entropic origin but it again leads to a reduction of late time diffusivity
[14, 15, 16, 17, 18]. The effect of the confining geometry on effective
diffusivity has been widely studied [19, 20, 2, 16, 21, 22, 23, 24, 25, 26,
27], but the vast majority of existing theories focuses on perfectly
reflecting boundaries.

However, the diffusion of a finite size tracer near a surface is actually much
more complicated that of simple Brownian motion at a reflecting wall. First,
the tracer will typically be subject to forces due to the surface, such as van
der Waals or electrostatic potentials [3]. Secondly, the components of the
diffusion tensor in the vicinity of the wall are reduced due to the no-slip
boundary condition on the ambient fluid [28]. A simple minimal model for these
effects is that of a tracer particle that can transiently attach (or adsorb)
to, diffuse on, and detach (desorb) from the surface, thus alternating between
phases of bulk diffusion or surface diffusion (Fig. 1). This type of model was
first proposed in the 90’s when it was realized that the effective lateral
dispersion of molecules adsorbed at a fluid-solid interface is significantly
modified by temporary excursions into the fluid [29]. The stochastic motion
resulting from interplay of bulk and surface diffusion has recently been
directly observed experimentally [30, 31, 32, 33, 34].

Figure 1: (color online) Examples of trajectories of a tracer particle
diffusing in the bulk (orange lines) or on the surface (thick blue lines) in a
regular array of spherical obstacles (a) or in a micro-channel (b).

It is also known that surface mediated transport gives rise to non-trivial
effects in the context of target search kinetics [35, 36, 37, 38] and that in
certain cases the search efficiency can be tuned via the parameters of the
surface mediated diffusion model. A classic related example is the
optimization of search on DNA by alternating phases of one dimensional
diffusion along the DNA, on which the target site is found, and three-
dimensional bulk excursions [39, 40, 41]. It is less clear if similar
optimization effects appear for transport properties. In Ref. [42], the effect
of adding a short range attraction to otherwise hard spheres was studied
numerically. Remarkably, it was found that spending time trapped at the
surface can lead to an increase in the late time diffusion of the tracer. The
study of [42] neglects the fact that diffusion near the hard sphere is slowed
down and one is naturally lead to ask the question as to whether the
enhancement effect persists when this is taken into account. Moreover, for
general geometries, there are no theories predicting this enhancement of
dispersion: existing approaches deal with uniform channels (where the effect
is absent) [43, 44] or in fast exchange limit for non-planar geometries [45,
4, 46], where no increase of dispersion due to sticky surfaces was found.

In this Letter, we introduce a theory that quantitatively predicts the
increase in late time dispersion (with respect to its value for reflecting
boundaries) induced by making the boundaries attractive, even when the surface
diffusivity is lower than the microscopic bulk diffusivity. As an example, we
give analytical formulas characterizing dispersion in non-dilute spherical
sticky obstacles, which explain the simulation results of [42]. Results are
then given for slowly undulating channels, in this case the increase of
dispersion turns out to be stronger. For strongly undulating channels, we show
that enhancement of dispersion no longer occurs. Our results identify in which
situations dispersion is optimized by the interplay between diffusion in bulk
and on the surface. The effect uncovered here can be viewed as an example of
catalysis for diffusion where the introduction of the surface state via
reaction with the surface induces an increase in the rate of dispersion.

Model of surface mediated transport - We consider a tracer particle, of
position $\mathbf{r}(t)$ at time $t$ in a $d$-dimensional space, an
heterogeneous medium with obstacles, or confining boundaries. The particle can
either diffuse in the bulk (“b”), with a bulk molecular diffusivity $D_{b}$,
and local drift $\mathbf{u}_{b}(\mathbf{r})$, or diffuse along the surface
(“s”) of the obstacles or confining surface, where it diffuses with diffusion
coefficient $D_{s}$ (typically smaller than $D_{b}$). We can also consider a
local drift field $\mathbf{u}_{s}(\mathbf{r})$, within the surface. By $k_{d}$
we denote a detachment rate at which the tracer desorbs from the surface. When
the transition between bulk and surface is viewed as a reaction, the binding
kinetics is quantified by an imperfect reactivity parameter $k_{a}$ [47],
which has the dimension of a velocity. In the context of Taylor dispersion, it
has been shown how the parameters $k_{a}$ and $k_{d}$ can be determined from a
microscopic model [48].

The probability densities $p_{b}(\mathbf{r},t)$ and $p_{s}(\mathbf{r},t)$ to
be at position $\mathbf{r}$ at time $t$ in the bulk and surface obey

$\displaystyle\frac{\partial p_{b}}{\partial t}$
$\displaystyle=\nabla\cdot[D_{b}\nabla
p_{b}-\textbf{u}_{b}(\mathbf{r})p_{b}],$ (1) $\displaystyle\frac{\partial
p_{s}}{\partial t}$
$\displaystyle=\nabla_{s}\cdot[D_{s}\nabla_{s}p_{s}-\textbf{u}_{s}(\mathbf{r})p_{s}]-k_{d}\
p_{s}+k_{a}\ p_{b},$ (2)

where $\nabla_{s}$ is the surface nabla operator. The boundary conditions,
determined from probability conservation, are given by

$\textbf{n}\cdot\left[D_{b}\nabla
p_{b}-\textbf{u}_{b}(\textbf{r})p_{b}\right]=k_{d}\ p_{s}-k_{a}\ p_{b},$ (3)

where n is the surface normal vector, pointing away from the bulk. Finally, we
assume that the geometry of the medium, as well as all fields
$\mathbf{u}_{b}(\mathbf{r}),\mathbf{u}_{s}(\mathbf{r})$, are spatially
periodic.

Effective transport properties -  The late time transport properties in a
given spatial direction (say the $x$ direction) are quantified by an effective
drift $v_{e}$ and effective diffusivity $D_{e}$ in this direction:

$\displaystyle v_{e}=\frac{\overline{x(t)-x(0)}}{t},\
D_{e}=\lim_{t\to\infty}\frac{\overline{[x(t)-x(0)-v_{e}t]^{2}}}{2t},$ (4)

where the overline denotes the average over stochastic trajectories. The
average drift is straightforward to calculate (see Appendix A)

$\displaystyle v_{e}=\int_{V}d\mathbf{r}\ (\nabla x)\cdot\mathbf{j}_{b}\
+\int_{S}dS\ (\nabla_{s}x)\cdot\mathbf{j}_{s},$ (5)
$\displaystyle\mathbf{j}_{b}=(\mathbf{u}_{b}-D_{b}\nabla)P_{b}^{st},\ \
\mathbf{j}_{s}=(\mathbf{u}_{s}-D_{s}\nabla_{s})P_{s}^{st},$ (6)

where the integrals are evaluated over the volume $V$ and the surface $S$ of
the walls in one period. The functions $P_{b}^{st}(\mathbf{r})$ and
$P_{s}^{st}(\mathbf{r})$ are the stationary probability densities of the
position modulo the period and are time-independent solutions of Eqs. (1) and
(2) with periodic conditions.

To compute $D_{e}$ in the absence of any advection, we use the fluctuation
dissipation theorem that relates it to the effective drift when a small force
$F$ acts on the tracer particle:

$\displaystyle D_{e}(F=0)=k_{B}T\times\left[\frac{d}{dF}\
v_{e}(F)\right]_{F=0},$ (7)

where $k_{B}T$ is the thermal energy. Here, a force $F$ in the $x$ direction
corresponds to the following drift fields:

$\displaystyle\mathbf{u}_{b}(\mathbf{r})=\beta D_{b}F(\nabla x),\
\mathbf{u}_{s}(\mathbf{r})=\beta D_{s}F(\nabla_{s}x),$ (8)

with $\beta=1/k_{B}T$. At zero force, we denote the stationary probability
densities by $P_{b}^{st,0},P_{s}^{st,0}$, both are uniform and are given by:

$\displaystyle P_{b}^{st,0}(\mathbf{r})=\frac{1}{V+S\delta},\
P_{s}^{st,0}(\mathbf{r})=\delta P_{b}^{st,0},$ (9)

where $\delta=k_{a}/k_{d}$ is an adsorption length (that can be much larger
than the range of interactions with the surface). We now introduce two
auxiliary fields $f_{b}$ and $f_{s}$ which quantify the deviation of the
stationary probability density from the uniform distribution at low forces:

$\displaystyle P_{b}^{st}\simeq P_{b}^{st,0}+\beta Ff_{b}(\mathbf{r}),\
P_{s}^{st}\simeq P_{s}^{st,0}+\beta Ff_{s}(\mathbf{r}).$ (10)

Inserting this ansatz into Eqs. (5) and (7), we see that $D_{e}(F=0)$ (denoted
simply by $D_{e}$) is given by

$\displaystyle D_{e}=$ $\displaystyle D_{b}\int_{V}\mbox{d}\textbf{r}\ (\nabla
x)\cdot\nabla\left[P_{b}^{st,0}x-f_{b}\right]$
$\displaystyle+D_{s}\int_{S}\mbox{d}S\
(\nabla_{s}x)\cdot\nabla_{s}\left[P_{s}^{st,0}x-f_{s}\right].$ (11)

The equations for $f_{b},f_{s}$ are obtained by inserting the ansatz (10) into
the transport equations (1, 2, 3) and expanding to linear order in $F$, which
yields

$\displaystyle\nabla^{2}f_{b}=0,\hskip
14.22636ptD_{s}\nabla_{s}^{2}[f_{s}-P_{s}^{st,0}x]=k_{d}f_{s}-k_{a}f_{b},$
(12) $\displaystyle D_{b}\ \textbf{n}\cdot\nabla
f_{b}=k_{d}f_{s}-k_{a}f_{b}+D_{b}P_{b}^{st,0}\ (\nabla x)\cdot\textbf{n}.$
(13)

The above equations, together with periodic boundary conditions, define
$f_{b}$ and $f_{s}$ up to an unimportant additive constant. This kind of
formula linking micro and macro-transport properties can also be found in
macrotransport theory [4], homogenization methods or Kubo formulas [49, 50].
In the case of surface mediated transport, the above equations are general and
in the limit of fast exchange between bulk and surface
($k_{d},k_{a}\to+\infty$ keeping $k_{a}/k_{d}$ constant) reproduce the results
given in [45, 4, 46]. Note that surface curvature effects arise in Eq. (12)
via the term $\nabla_{s}^{2}x$ which does not vanish (as it would in the
bulk). The above partial differential equations can be numerically solved by
standard finite element routines. However, we will also derive analytical
results in certain limits.

Figure 2: (color online) Long-time dispersion for a periodic array of
spherical obstacles ($d=3$) with $D_{s}=D_{b}$. Lines: theoretical prediction,
Eq. (15). Squares: data obtained by numerical integration of Eqs. (12), (13).
Triangles: results of stochastic simulations of the stochastic process
$\mathbf{r}(t)$ associated to the Fokker-Planck equations (1,2). Circles:
Monte Carlo simulation data taken from [42], where the tracer particles are
submitted to an exponentially decaying potential. The parameter $\delta$ was
adjusted so that the stationary probability densities (9) in our description
are the same as in the case of this potential (see Appendix D.3.a. for
details).

Dispersion in regular arrays of spheres -  One the most studied models of a
crowded environment is that of a regular array of spherical obstacles. Here we
consider the cases of a square lattice (in 2D, $d=2$) or a cubic lattice (in
3D, $d=3$), see Fig. 1(a). We denote by $L$ the distance between the centers
of nearest neighbouring spheres, $R$ their radius, and $\varphi$ their volume
fraction. An exact analytical expression for the dispersion for any $\varphi$
is difficult to obtain even without the surface interactions. However, in the
limit of small $R/L$ (equivalently, the small $\varphi$ limit), we can look
for solutions in terms of matched asymptotic expansions, with an inner
solution for $f_{b}(\mathbf{r})$ that varies at the scale $|\mathbf{r}|\sim R$
, and an outer solution varying at the scale $L$. The auxiliary fields in the
limit $R\to 0$ are written as

$\displaystyle f_{b}(\mathbf{r})=\begin{cases}\sum\limits_{n\geq
0}R^{n+1}f_{b}^{(n)}(r/R,\omega)&\left(r\ll L,\text{inner}\right)\\\
\sum\limits_{n\geq 0}R^{d+n}F_{b}^{(n)}(r,\omega)&(r\gg
R,\text{outer})\end{cases},$ (14)

where $r$ is the distance to the center of the nearest obstacle and $\omega$
denote angular variables. The solutions at successive orders are found by
inserting Eq. (14) into our formalism and imposing that both outer and inner
solutions are equal in the regime $R\ll r\ll L$ (where they are both valid).
To order $\varphi^{4}$ (see Appendix D), we find

$\displaystyle D_{e}=$
$\displaystyle\frac{D_{b}}{1+(d\xi-1)\varphi}\frac{1-\alpha\varphi}{1+\alpha\varphi/(d-1)}+o(\varphi^{4}),$
$\displaystyle\alpha=$
$\displaystyle\frac{D_{b}+D_{s}\xi(\gamma-1)(d-1)}{D_{b}+D_{s}\xi[\gamma(d-1)+1]},\
\gamma=\frac{D_{b}}{k_{a}R},\ \xi=\frac{\delta}{R}.$ (15)

The above formula is exact up to the order $\varphi^{4}$ in the small
$\varphi$ limit, when $\xi$ and $\gamma$ are kept constant. It is valid for a
large range of values of $\varphi$, as shown in Fig. 2 where it is compared to
numerical solution of the full transport equations (12, 13). The theory also
agrees with the Monte Carlo simulations of Ref. [42]. We can define a critical
surface diffusivity $D_{s}^{*}$ as the value of $D_{s}$ above which a weak
attraction increases $D_{e}$ (so $\left[\partial
D_{e}/\partial\delta\right]_{\delta=0}=0$ at $D_{s}=D_{s}^{*}$). Using Eq.
(15), we find

$\displaystyle D_{s}^{*}=D_{b}\left(1-\frac{1-\varphi}{d}\right).$ (16)

Interestingly, $D_{s}^{*}$ depends only on $\varphi$ and is independent of the
adsorption and desorption rates. At this level of perturbation theory we also
see that $D_{s}^{*}<D_{b}$. Consequently, if the surface and bulk
diffusivities are equal ($D_{s}=D_{b}$), making the surfaces weakly attractive
generically increases dispersion. This effect holds even when $D_{s}<D_{b}$
(as long as $D_{s}>D_{s}^{*}$). Hence, spending time on a surface, even with
reduced diffusion, can actually lead to faster dispersion than in the case of
purely reflecting obstacles. A similar formula to Eq. (15) was found in Ref.
[51] in the fast exchange limit ($\gamma=0$), however this formula does not
predict an increase of dispersion for weakly attractive spheres, and it
disagrees with our theory and with numerical simulations (see Appendix
D.3.b.). Finally, the increase in $D_{e}$ is rather marginal in this geometry,
as can be seen in Fig. 2.

Dispersion in slowly undulating channels -  We now analyze the case of a
surface-mediated diffusion process in a symmetric channel of half-width $h(x)$
(if $d=2$) or in an axisymmetric channel of radius $h(x)$ (if $d=3$), see
figure 1(b). Such geometries are often viewed as a paradigm for transport in
confined media [14, 17]. We consider the limit of slowly varying undulations
($L\to\infty$), where the problem can be solved using a perturbation
expansion:

$\displaystyle f_{b}=\sum_{n\geq
0}\frac{1}{L^{n}}f_{b}^{(n)}\left(\frac{x}{L},\mathbf{r}_{\perp}\right),f_{s}=\sum_{n\geq
0}\frac{1}{L^{n}}f_{s}^{(n)}\left(\frac{x}{L}\right),$ (17)

where $\mathbf{r}_{\perp}$ is the position perpendicular to the direction $x$.
Inserting the above ansatz into our formalism and writing all equations order
by order, we obtain (see Appendix E)

$\displaystyle D_{e}\underset{L\to\infty}{=}\frac{1}{\left\langle
h^{d-2}\left(h+(d-1)\delta\right)\right\rangle\left\langle\frac{h^{2-d}}{D_{b}h+(d-1)D_{s}\delta}\right\rangle},$
(18)

where $\left<\cdot\right>=\frac{1}{L}\int_{0}^{L}\cdot\ dx$ denotes the
spatial average over one period $L$. Note that $D_{e}$ is not given by a
simple steady state average between effective bulk diffusion and surface
diffusion as it is the case for a flat channel [43]. For reflecting walls
($\delta=0$), we recover the well known result obtained with the Fick-Jacobs
approximation, $D_{e}(\delta=0)=D_{b}[\langle h^{d-1}\rangle\langle
h^{-(d-1)}\rangle]^{-1}$ [16, 24]. Once again, defining $D_{s}^{*}$ as the
value of $D_{s}$ for which $[\partial D_{e}/\partial\delta]_{\delta=0}=0$, we
find

$\displaystyle D_{s}^{*}=D_{b}\frac{\langle h^{d-2}\rangle\langle
h^{1-d}\rangle}{\langle h^{d-1}\rangle\langle h^{-d}\rangle}.$ (19)

It can be shown that $D_{s}^{*}$ is smaller than $D_{b}$ for any choice of
profile $h$. Figure 3 illustrates the increase of dispersion when surfaces are
made weakly attractive for one example of a two-dimensional channel.
Interestingly, the increase of dispersion can be as high as $250\%$ (upper
curve), this change could be made arbitrarily high for sharper channel
undulations. The enhancement effect can thus be significantly larger than that
observed for sticky spheres as studied above and in Ref. [42].

It is clear that the enhancement of diffusion is related to the fact that
diffusing along the surface avoids bottlenecks (or entropic barriers) in the
channel. However, the effect of the surface interaction is more subtle than
just this basic effect. As seen in Figure 3 from the curve with solid circles,
a finite interaction with the surface can enhance diffusion with respect to
both the case $\delta=0$ (pure reflection) and $\delta=\infty$ (diffusion
restricted to the surface), even when the surface diffusion is $30\%$ of that
of bulk diffusion. Actually, it can be shown that when
$D_{e}(0)=D_{e}(\infty)$, then the inequality $D_{s}>D_{s}^{*}$ always holds
(see Appendix E.3). Hence, for any channel geometry, there is a regime for
which the effective diffusivity is larger than the effective diffusivity for
completely reflecting ($\delta=0$) or completely sticky ($\delta=\infty$)
boundaries. The subtlety of the effect is also highlighted by the earlier
results on spherical obstacles, where the enhancement effect actually shows up
at first order in the volume fraction, thus in absence of significant
bottlenecks.

Figure 3: (color online) Long-time diffusion coefficient for a two-dimensional
channel of height profile $h/a=16+15\arctan[\cos u/(\sin
u+3/2)]/\arctan(2/\sqrt{5})$ with $u=2\pi x/L$. Symbols: results of the
numerical integration of Eqs. (12) and (13) for $k_{a}a/D_{b}=10$ and
$L/a=500$ ; for this value the channel shape is shown in the inset. Lines:
predictions of Eq. (18) in the limit of slowly varying width.

Dispersion in highly undulating channels - Finally, we investigate rapidly
varying channels ($L\ll a$, with $a$ the minimal channel height), for which we
obtain (see Appendix F)

$\displaystyle D_{e}$ $\displaystyle\underset{L\ll
a}{=}\frac{V_{c}}{V+S\delta}D_{b},$ (20)

where $V_{c}$ is the volume of a central region (a cylinder of radius $a$, so
$V_{c}=2aL$ in 2D and $V_{c}=\pi a^{2}L$ in 3D). This formula can be deduced
from a simple ergodic argument. In the central region, the tracer particle can
freely diffuse, but in the peripheral regions made of very thin dead-ends, the
tracer particle does not contribute to dispersion along the channel axis (at
leading order). Hence, the late time diffusivity is the product of $D_{b}$
with the fraction of time spent in the central region, leading to Eq. (20). In
this geometric limit, we see that surface interactions only decrease
dispersion.

Conclusion -  It is well established that first passage kinetics to a target
can be optimized by an appropriate combination of sejourns of diffusion in
spaces of different dimensions [36, 40, 35], as exemplified by the search of a
gene sequence on DNA [39, 41, 40]. However, this optimization is not linked to
an acceleration of diffusivity. Here, we have theoretically demonstrated a
similar optimization effect holds for transport features: dispersion in
crowded media can be enhanced if the obstacles or surfaces exhibit a short
ranged attraction for the tracer. This enhancement phenomenon was first noted
in Ref. [42] for hard spheres. Using a surface mediated diffusion model, we
have developed a general transport theory to analyze the effect of short range
attractive interactions on dispersion. This theory explains the results of
Ref. [42] and we have provided analytical results for dilute systems (which
work rather well even for large volume fractions). We have also analyzed
dispersion in symmetric channels where the particle can adsorb and detach from
the walls. For slowly undulating channels, we have shown that an even stronger
diffusion enhancement can occur and have given exact results. In both
geometries, we find that even when the surface diffusion constant $D_{s}$ is
smaller than that in the bulk $D_{b}$ (as must be the case physically) this
enhancement can still occur and we have determined the critical values of
$D_{s}$ above which this happens.

##  Appendix A Calculation of the effective drift - derivation of Eq. (5)

Here, we compute the effective drift defined in Eq. (4). For now, we consider
a fixed initial position $\mathbf{r}_{0}$ and initial state $i_{0}$ (equal to
$b$ or $s$, i.e. in the bulk or on the surface). The displacement in the $x$
direction after a time $t$ reads

$\overline{x(t)-x_{0}}=\int_{\Omega}d\mathbf{r}(x-x_{0})p_{b}(\mathbf{r},t|\mathbf{r}_{0},i_{0})+\int_{\Sigma}dS(\mathbf{r})(x-x_{0})p_{s}(\mathbf{r},t|\mathbf{r}_{0},i_{0}),$
(21)

where $\Omega$ and $\Sigma$ are the volume and surface of the medium. Taking
the temporal derivative and using Eqs. (1,2) yields

$\displaystyle\partial_{t}\overline{x(t)}=$
$\displaystyle\int_{\Omega}d\mathbf{r}\
(x-x_{0})\nabla\cdot\left[(D_{b}\nabla-\mathbf{u}_{b})p_{b}(\mathbf{r},t|\mathbf{r}_{0},i_{0})\right]$
$\displaystyle+\int_{\Sigma}dS(\mathbf{r})(x-x_{0})\left\\{\nabla_{s}\cdot\left[(D_{s}\nabla_{s}-\mathbf{u}_{s})p_{s}(\mathbf{r},t|\mathbf{r}_{0},i_{0})\right]-k_{d}\
p_{s}(\mathbf{r},t|\mathbf{r}_{0},i_{0})+k_{a}\
p_{b}(\mathbf{r},t|\mathbf{r}_{0},i_{0})\right\\}.$ (22)

After performing integrations by parts and using the boundary condition Eq.
(3), we obtain

$\displaystyle\partial_{t}\overline{x(t)}=\int_{\Omega}d\mathbf{r}(\nabla
x)\cdot\left(\mathbf{u}_{b}-D_{b}\nabla\right)p_{b}(\mathbf{r},t|\mathbf{r}_{0},i_{0})+\int_{\Sigma}dS(\mathbf{r})(\nabla_{s}x)\cdot\left(\mathbf{u}_{s}-D_{s}\nabla_{s}\right)p_{s}(\mathbf{r},t|\mathbf{r}_{0},i_{0}).$
(23)

We now introduce $P_{s}(\mathbf{r},t)$ and $P_{b}(\mathbf{r},t)$, which are
the probability densities at $t$ to be respectively on the surface or in the
bulk, at position $\mathbf{r}$ modulo the period, i.e.:

$\displaystyle
P_{b}(\mathbf{r},t)=\sum_{\mathbf{L}}p_{b}(\mathbf{r}+\mathbf{L},t),\
P_{s}(\mathbf{r},t)=\sum_{\mathbf{L}}p_{s}(\mathbf{r}+\mathbf{L},t),$ (24)

where the vectors $\mathbf{L}$ are all the vectors such that the environments
at $\mathbf{r}$ and $\mathbf{r}+\mathbf{L}$ are identical (for example, in the
channel geometry $\mathbf{L}=nL\mathbf{e}_{x}$ with $n$ integer). In Eq. (23),
$\nabla x$, $\nabla_{s}x$, $\mathbf{u}_{b}$ and $\mathbf{u}_{s}$ are periodic
in $\mathbf{r}$. Therefore, splitting the integral into multiple integrals
(each of them covering one period) and using the above relation yields

$\displaystyle\partial_{t}\overline{x(t)}=\int_{V}d\mathbf{r}\ (\nabla
x)\cdot\left(\mathbf{u}_{b}-D_{b}\nabla\right)P_{b}(\mathbf{r},t|\mathbf{r}_{0},i_{0})+\int_{S}dS(\mathbf{r})\
(\nabla_{s}x)\cdot\left(\mathbf{u}_{s}-D_{s}\nabla_{s}\right)P_{s}(\mathbf{r},t|\mathbf{r}_{0},i_{0}),$
(25)

where $V$ and $S$ are the volume and surface in one period of the channel
only. This result is obtained by assuming a fixed initial state and position,
let us now extend it to the case that $\mathbf{r}_{0}$ and $i_{0}$ are drawn
from the stationary distribution of the dynamics modulo the period. In this
case, it is useful to note that

$P_{i}^{st}(\mathbf{r})=\int_{V}d\mathbf{r}_{0}P_{i}(\mathbf{r},t|\mathbf{r}_{0},b)P_{b}^{st}(\mathbf{r}_{0})+\int_{S}dS(\mathbf{r}_{0})P_{i}(\mathbf{r},t|\mathbf{r}_{0},s)P_{s}^{st}(\mathbf{r}_{0})\
\ ,\forall t\hskip 56.9055pt(i\in\\{b,s\\}),$ (26)

so that averaging Eq. (25) over initial conditions and integrating over $t$
leads to

$\displaystyle
v_{e}=\partial_{t}\overline{x(t)}=\frac{\overline{x(t)-x(0)}}{t}=\int_{V}d\mathbf{r}\
(\nabla
x)\cdot\left(\mathbf{u}_{b}-D_{b}\nabla\right)P_{b}^{st}(\mathbf{r})+\int_{S}dS(\mathbf{r})\
(\nabla_{s}x)\cdot\left(\mathbf{u}_{s}-D_{s}\nabla_{s}\right)P_{s}^{st}(\mathbf{r}),$
(27)

which is Eq. (5) in the main text. Note that this expression also holds, at
long times only, for any choice of initial condition, as can be seen by taking
the long time limit of Eq. (25).

##  Appendix B Tensor analysis for surface diffusion

We recall some basic results of tensor analysis that are useful to interpret
explicitly $\nabla_{s}^{2}x$ in different coordinate systems. We assume that
the surface in a $d$ dimensional space is parametrized by the coordinates
$y^{i}$ ($1\leq i\leq d-1$). If $\mathbf{r}_{s}$ is the position vector on the
surface, then the metric given by

$g_{ij}=\textbf{e}_{s,i}\cdot\textbf{e}_{s,j},\hskip
28.45274pt\textbf{e}_{s,i}=\frac{\partial\textbf{r}_{s}}{\partial y^{i}},$
(28)

where the suffix $s$ indicates that the position remains on the surface. We
use the notation $|g|=\det(g_{ij})$ and $g^{ij}$ indicates the inverse of the
metric tensor (i.e. $g^{ij}g_{jk}=\delta_{ik}$, with Einstein’s summation
convention). In our formalism, we need to identify the elementary surface
elements, the surface gradient and Laplacian of a scalar function $\phi$,
which are well known and are given by

$dS=dy^{1}dy^{2}...dy^{d-1}\sqrt{|g|},\hskip 28.45274pt\nabla_{s}\phi=g^{ij}\
\frac{\partial\phi}{\partial y^{j}}\ \textbf{e}_{s,i},\hskip
28.45274pt\nabla_{s}^{2}\phi=\frac{1}{\sqrt{|g|}}\frac{\partial}{\partial
y^{i}}\left(\sqrt{|g|}\ g^{ij}\ \frac{\partial\phi}{\partial y^{j}}\right).$
(29)

##  Appendix C Algorithm for stochastic simulations

We have used the following algorithm to simulate the stochastic process
$\mathbf{r}(t)$ associated to the Fokker-Planck equations (1),(2) in periodic
arrays of spheres (the algorithm could be adapted to other geometries):

  * •

a particle in the bulk evolves with $dr_{i}=\sqrt{2D_{b}}dW_{i}$ (with
$\langle dW_{i}\rangle=0$ and $\langle dW_{i}dW_{j}\rangle=dt\delta_{ij}$)
during the time step $dt$ in each direction $i=\\{x,y,z\\}$.

  * •

At the end of a step, if the tracer particle end up inside the obstacle, it
switches to the attached state with probability
$P_{\text{adsorb}}=k_{a}\sqrt{\pi dt/D_{b}}$ and it is reflected with
probability $1-P_{\text{adsorb}}$, see Ref. [53].

  * •

If the particle is attached to a spherical 3D obstacle, the stochastic
equation (with Ito convention) for the angular coordinates $(\theta,\phi)$ is

$d\theta_{t}=\frac{D_{s}dt}{R^{2}\tan\theta_{t}}+\frac{\sqrt{2D_{s}}}{R}dW_{\theta}\
,\ d\phi_{t}=\sqrt{\frac{2D_{s}}{R^{2}\sin^{2}\theta_{t}}}dW_{\phi},\
dW_{\phi}^{2}=dW_{\theta}^{2}=dt.$ (30)

  * •

An attached particle switches to the detached state with probability
$1-e^{-k_{d}\cdot dt}$.

##  Appendix D Dispersion in periodic array of spheres - derivation of Eq.
(15)

###  D.1 2D obstacles

Here, we compute $D_{e}$ for a square array of 2D spheres of radius $R$. We
use polar coordinates $(r,\theta)$, with $r$ the distance to the center of the
disk and $\theta$ the angle with the $x$-axis. Here we chose the units of
length so that $L=1$, without loss of generality. The surface is parametrized
by the coordinate $y^{1}=\theta$. Using the formulas of Appendix B, we obtain
$\nabla_{s}^{2}x=-\cos\theta/R$, and the equations in polar coordinates are:

$\left\\{\begin{split}&\text{(bulk) }\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial f_{b}}{\partial
r}\right)+\frac{1}{r^{2}}\frac{\partial^{2}f_{b}}{\partial\theta^{2}}=0,&(r>R)\\\
&\text{(surface)
}\frac{D_{s}}{R^{2}}\frac{\partial^{2}f_{s}}{\partial\theta^{2}}=k_{d}f_{s}(\theta)-k_{a}f_{b}(r=R,\theta)-D_{s}P_{s}^{st,0}\frac{\cos\theta}{R},&(r=R)\\\
&\text{(boundary condition) }-D_{b}\frac{\partial f_{b}}{\partial
r}(r=R,\theta)=k_{d}f_{s}-k_{a}f_{b}-D_{b}P_{b}^{st,0}\cos\theta.&(r=R)\end{split}\right.$
(31)

After an integration by parts in Eq. (11), we find that the effective
diffusion constant is

$D_{e}=D_{b}(1-\varphi)P_{b}^{st,0}+D_{s}\varphi\xi
P_{b}^{st,0}+D_{b}R\int_{0}^{2\pi}\mbox{d}\theta\
\cos(\theta)f_{b}(r=R,\theta)-D_{s}\int_{0}^{2\pi}\mbox{d}\theta\cos(\theta)f_{s}(\theta).$
(32)

We introduce the dimensionless parameters $\xi=\delta/R$ and
$\gamma=D_{b}/(k_{a}R)$. We now compute $D_{e}$ in the limit $R\to 0$ (or,
equivalently, the limit of vanishing volume fraction $\varphi=\pi R^{2}\to
0$), when $\xi$ and $\gamma$ are kept constant. We solve this problem using
strongly localized perturbation analysis [54], by considering that the
solution $f_{b}$ varies at two different length scales $r\sim R$ and $r\sim 1$
when $R\to 0$ (we recall that the period is $L=1$). We introduce an inner
solution (valid near the obstacles), depending on $\tilde{r}=r/R$ and
$\theta$, which has an expansion in the limit $R\rightarrow 0$ at fixed
$\tilde{r}=r/R$:

$\begin{split}f_{b}(\mathbf{r})=P_{b}^{st,0}R\left[f_{b}^{0}(\tilde{r},\theta)+R^{2}f_{b}^{1}(\tilde{r},\theta)+R^{4}f_{b}^{2}(\tilde{r},\theta)+...\
\right],f_{s}(\mathbf{r})=P_{b}^{st,0}R^{2}\left[f_{s}^{0}(\theta)+R^{2}f_{s}^{1}(\theta)+R^{4}f_{s}^{2}(\theta)+...\
\right].\end{split}$ (33)

Note that we have introduced the factor $P_{b}^{st,0}=1/(1+(2\xi-1)\varphi)$
since it is obvious in Eq. (32) that the solutions are proportional to this
factor. We introduce an outer solution, far from the obstacle, which admits
the following expansion when $R\rightarrow 0$ at fixed $r$:

$f_{b}(\mathbf{r})=P_{b}^{st,0}R^{2}\left[F_{b}^{0}(r,\theta)+R^{2}F_{b}^{1}(r,\theta)+R^{4}F_{b}^{2}(r,\theta)+...\
\right].$ (34)

It follows that all $f_{b}^{i}(\tilde{r},\theta)$ and $F_{b}^{i}(r,\theta)$
are harmonic functions for all orders $i$, and that all $F_{b}^{i}(r,\theta)$
must satisfy periodic conditions on the sides of the unit square. The
equations for $f_{s}^{i}$ and for the value of $f_{b}^{i}$ near the obstacles
are found by inserting Eq. (33) into Eq. (31) and identifying terms order by
order, leading to:

$D_{s}\gamma\xi\frac{\partial^{2}f_{s}^{i}}{\partial\theta^{2}}=D_{b}\left(f_{s}^{i}-\xi
f_{b}^{i}\right)-\delta_{i,0}\xi^{2}D_{s}\gamma\cos\theta,\hskip
28.45274pt-\xi\gamma\frac{\partial
f_{b}^{i}}{\partial\tilde{r}}(\tilde{r}=1,\theta)=f_{s}^{i}-\xi
f_{b}^{i}-\delta_{i,0}\xi\gamma\cos\theta.$ (35)

Additional conditions to determine $f_{b}^{i}$ and $F_{b}^{i}$ come from the
requirement that the inner and the outer solution coincide in their regime of
common validity $R\ll r\ll 1$, so that the expansion of Eq. (33) for
$\tilde{r}\to\infty$ coincides with the expansion of (34) for $r\to 0$.

At leading order, $f_{b}^{0}$ is the harmonic function that does not diverge
for large $\tilde{r}$ and that satisfies Eq. (35) for $i=0$:

$\displaystyle
f_{b}^{0}(\tilde{r},\theta)=\underbrace{\frac{A_{0}}{\tilde{r}}\cos(\theta)}_{\text{must
match }F_{b}^{0}},\hskip
8.5359ptf_{s}^{0}(\theta)=\tilde{f}_{s}^{0}\cos(\theta),\hskip
8.5359ptA_{0}=-\alpha,\hskip 8.5359pt\tilde{f}_{s}^{0}=\xi\
\frac{D_{s}\xi(\gamma+1)-D_{b}}{D_{s}\xi(\gamma+1)+D_{b}},\hskip
8.5359pt\alpha=\frac{D_{b}+D_{s}\xi(\gamma-1)}{D_{b}+D_{s}\xi(\gamma+1)}.$
(36)

We see that $f_{b}^{0}\approx A_{0}\cos\theta/\tilde{r}$ for large
$\tilde{r}$, so that $f_{b}\simeq A_{0}\cos\theta R^{2}/r$ in the region $R\ll
r\ll 1$. This implies that the leading order of the outer solution must be of
order $R^{2}$, as assumed in Eq. (34), and that $F_{b}^{0}(r\to
0,\theta)\simeq A_{0}\cos\theta/r$, as indicated in Eq. (36). Since
$F_{b}^{0}$ is also harmonic with periodic boundary conditions, it can be
expressed as

$\displaystyle F_{b}^{0}(r,\theta)=-2\pi A_{0}\ \textbf{e}_{x}\cdot\nabla G,$
(37)

where $G$ is the pseudo Green’s function for the unit-square with periodic
conditions, which satisfies the equation:

$-\nabla^{2}G=\delta(x)\delta(y)-1,\ G(x\pm 1,y)=G(x,y\pm 1)=G(x,y).$ (38)

The Green’s function $G$ is given in closed form in Ref. [55], its behaviour
near the origin reads

$G(r,\theta)=-\frac{\log
r}{2\pi}+\frac{r^{2}}{4}+C-\frac{C_{0}}{8\pi}r^{4}\cos(4\theta)+\mathcal{O}(r^{8}),$
(39)

with $C_{0}$ a constant that we will not need to detail here. Eq. (37) is
justified by noting that $F_{b}^{0}$ near the origin is

$\displaystyle
F_{b}^{0}(r,\theta)=\underbrace{\frac{A_{0}}{r}\cos(\theta)}_{\text{ matching
}f_{b}^{0}}\underbrace{-\pi A_{0}\cos(\theta)r}_{\text{must match
}f_{b}^{1}}+\underbrace{C_{0}\cos(3\theta)r^{3}}_{\text{must match
}f_{b}^{2}}+\mathcal{O}(r^{7}),$ (40)

where we see that the $1/r$ term correctly matches the corresponding
$1/\tilde{r}$ term in Eq. (36). We have also indicated that the terms of order
$r$ and $r^{3}$ will have to be matched by corresponding terms $\tilde{r}$ and
$\tilde{r}^{3}$ in the expansions of $f_{b}^{1}$ and $f_{b}^{2}$,
respectively. In particular, $f_{b}^{1}(\tilde{r},\theta)\simeq-\pi
A_{0}\tilde{r}\cos\theta$ for large $\tilde{r}$, since $f_{b}^{1}$ is
harmonic, it is of the form

$\displaystyle
f_{b}^{1}(\tilde{r},\theta)=\underbrace{\frac{A_{1}}{\tilde{r}}\cos(\theta)}_{\text{must
match }F_{b}^{1}}+\underbrace{B_{1}\tilde{r}\cos(\theta)}_{\text{matching
}F_{b}^{0}},\hskip 28.45274ptB_{1}=-\pi A_{0},$ (41)

Solving Eq. (35) for $i=1$, we obtain $f_{s}^{1}$ and the constant $A_{1}$:

$\displaystyle A_{1}=\alpha B_{1},\hskip
28.45274ptf_{s}^{1}(\theta)=\tilde{f}_{s}^{10}\cos(\theta),\hskip
28.45274pt\tilde{f}_{s}^{10}=\frac{2\xi
D_{b}B_{1}}{D_{s}\xi(\gamma+1)+D_{b}}.$ (42)

The outer solution at next-to-leading order $F_{b}^{1}$ and its expansion near
the origin then reads:

$F_{b}^{1}(r,\theta)=-2\pi A_{1}\textbf{e}_{x}\cdot\nabla
G=\underbrace{\frac{A_{1}}{r}\cos(\theta)}_{\text{ matching
}f_{b}^{1}}\underbrace{-\pi A_{1}\cos(\theta)r}_{\text{must match
}f_{b}^{2}}+\underbrace{C_{1}\cos(3\theta)r^{3}}_{\text{must match
}f_{b}^{3}}+\mathcal{O}(r^{7}).$ (43)

Second order: We continue the iteration and we look for $f_{b}^{2}$ in the
form

$f_{b}^{2}(\tilde{r},\theta)=\underbrace{\frac{A_{2}}{\tilde{r}}\cos(\theta)}_{\text{must
match }F_{b}^{2}}+\underbrace{B_{2}\tilde{r}\cos(\theta)}_{\text{matching
}F_{b}^{1}}+\underbrace{\frac{K_{2}}{\tilde{r}^{3}}\cos(3\theta)}_{\text{must
match }F_{b}^{3}}+\underbrace{C_{0}\tilde{r}^{3}\cos(3\theta)}_{\text{matching
}F_{b}^{0}},\hskip 28.45274ptB_{2}=-\pi A_{1},$ (44)

so it matches corresponding terms for $F_{b}^{1}$ and $F_{b}^{0}$ in Eqs. (40)
and (43), respectively. The expression for $A_{2}$ and ${f}_{s}^{2}$ can be
found by solving Eq. (35) with $i=2$:

$A_{2}=\alpha B_{2},\
f_{s}^{2}(\theta)=\tilde{f}_{s}^{20}\cos(\theta)+\tilde{f}_{s}^{21}\cos(3\theta),\
\tilde{f}_{s}^{20}=\frac{2\xi D_{b}B_{2}}{D_{s}\xi(\gamma+1)+D_{b}}.$ (45)

The values of $K_{2}$ and $\tilde{f}_{s}^{21}$ are not written because they do
not contribute to $D_{e}$. For the outer solution, we obtain

$F_{b}^{2}(r,\theta)=-2\pi A_{2}\textbf{e}_{x}\cdot\nabla
G=\underbrace{\frac{A_{2}}{r}\cos(\theta)}_{\text{ matching
}f_{b}^{2}}\underbrace{-\pi A_{2}\cos(\theta)r}_{\text{must match
}f_{b}^{3}}+\underbrace{C_{2}\cos(3\theta)r^{3}}_{\text{must match
}f_{b}^{4}}+\mathcal{O}(r^{7}).$ (46)

Third order: At this order, the expression of the inner solution reads

$f_{b}^{3}(\tilde{r},\theta)=\underbrace{\frac{A_{3}}{\tilde{r}}\cos(\theta)}_{\text{must
match }F_{b}^{3}}+\underbrace{B_{3}\tilde{r}\cos(\theta)}_{\text{matching
}F_{b}^{2}}+\underbrace{\frac{K_{3}}{\tilde{r}^{3}}\cos(3\theta)}_{\text{must
match }F_{b}^{4}}+\underbrace{C_{1}\tilde{r}^{3}\cos(3\theta)}_{\text{matching
}F_{b}^{1}},\hskip 28.45274ptB_{3}=-\pi A_{2}.$ (47)

This form is chosen so that the corresponding terms in $F_{b}^{2}$ and
$F_{b}^{1}$ are matched. Identically to Eq. (45), we get

$A_{3}=\alpha B_{3},\
f_{s}^{3}(\theta)=\tilde{f}_{s}^{30}\cos(\theta)+\tilde{f}_{s}^{31}\cos(3\theta),\
\tilde{f}_{s}^{30}=\frac{2\xi D_{b}B_{3}}{D_{s}\xi(\gamma+1)+D_{b}}.$ (48)

At this order of perturbation, we find that the effective diffusivity as given
in Eq. (32) is given by:

$\displaystyle
D_{e}=P_{b}^{st,0}\left\\{D_{b}(1-\varphi)+D_{s}\xi\varphi+D_{b}\varphi\left[A_{0}+\sum_{i=1}^{3}R^{2i}(A_{i}+B_{i})\right]-D_{s}\varphi\sum_{i=0}^{3}R^{2i}\tilde{f}_{s}^{i}\right\\}+o(R^{8}).$
(49)

An explicit calculation of all terms, followed by a re-summation, leads to the
result for $D_{e}$ for 2D circular obstacles:

$D_{e}=D_{b}P_{b}^{st,0}\left(1-2\alpha\varphi+2\alpha^{2}\varphi^{2}-2\alpha^{3}\varphi^{3}+2\alpha^{4}\varphi^{4}\right)+o(\varphi^{4})=\frac{D_{b}}{1+(2\xi-1)\varphi}\
\frac{1-\alpha\varphi}{1+\alpha\varphi}+o(\varphi^{4}).$ (50)

###  D.2 3D spherical obstacles

We now consider a 3-dimensional cubic array of spheres of radius $R$. Here, we
use spherical coordinates $(r,\theta,\phi)$, and we calculate the late time
diffusion constant in the $z$ direction, the distance between neighboring
spheres is again taken as $L=1$. Note that $r=0$ corresponds to the center of
an obstacle, and the nearest obstacles point to the directions
$(\theta=\pi/2,\phi=n\pi/2$), and $\theta=0$ and $\theta=\pi$. The long-time
diffusion constant is given by

$D_{e}=P_{b}^{st,0}\left[D_{b}\left(1-\varphi\right)+2\varphi\xi
D_{s}\right]+\int_{0}^{\pi}\mbox{d}\theta\int_{0}^{2\pi}\mbox{d}\phi\sin(\theta)\cos(\theta)\left(D_{b}R^{2}f_{b}-2D_{s}Rf_{s}\right),$
(51)

where $\varphi=4\pi R^{3}/3$ is the volume fraction. We adapt the method
exposed in the 2D case, by looking for an inner solution, for which the ansatz
reads:

$\left\\{\begin{split}f_{b}(\mathbf{r})&=P_{b}^{st,0}R\left[f_{b}^{0}(\tilde{r},\theta,\phi)+Rf_{b}^{1}(\tilde{r},\theta,\phi)+R^{2}f_{b}^{2}(\tilde{r},\theta,\phi)+...\
\right],\\\
f_{s}(\mathbf{r})&=P_{b}^{st,0}R^{2}\left[f_{s}^{0}(\theta,\phi)+Rf_{s}^{1}(\theta,\phi)+R^{2}f_{s}^{2}(\theta,\phi)+...\
\right],\end{split}\right.\hskip 28.45274pt(R\to 0,\ \text{fixed }\
\tilde{r}=r/L,\theta,\phi)$ (52)

where we have factorized by
$P_{b}^{st,0}=\left[1+(3\xi-1)\varphi\right]^{-1}$. The outer solution is

$f_{b}(\mathbf{r})=P_{b}^{st,0}R^{3}\left[F_{b}^{0}(r,\theta,\phi)+RF_{b}^{1}(r,\theta,\phi)+R^{2}F_{b}^{2}(r,\theta,\phi)+...\
\right].\hskip 28.45274pt(R\to 0,\ \text{fixed }\ r,\theta,\phi)$ (53)

Note that we have included all powers of $R$ in our ansatz although many of
these terms will be found to be zero afterwards. The solutions at successive
orders are found by requiring that (i) all $f_{b}^{i}$ and $F_{b}^{i}$ are
harmonic functions, (ii) $F_{b}^{i}$ satisfy periodic boundary conditions,
(iii) the behavior of the inner and the outer solution is the same at all
orders in the region $R\ll r\ll 1$, and (iv) the inner functions satisfy

$\left\\{\begin{split}&\text{(surface) }\gamma
D_{s}\xi\left[\frac{1}{\sin(\theta)}\frac{\partial}{\partial\theta}\left(\sin(\theta)\frac{\partial
f_{s}^{i}}{\partial\theta}\right)+\frac{1}{\sin(\theta)^{2}}\frac{\partial^{2}f_{s}^{i}}{\partial\phi^{2}}\right]=D_{b}\left(f_{s}^{i}-\xi
f_{b}^{i}\right)-2\delta_{i,0}\xi^{2}D_{s}\gamma\cos(\theta),\\\
&\text{(boundary condition) }-\xi\gamma\frac{\partial
f_{b}^{i}}{\partial\tilde{r}}(\tilde{r}=1,\theta)=f_{s}^{i}-\xi
f_{b}^{i}-\delta_{i,0}\xi\cos(\theta).\end{split}\right.$ (54)

At leading order, the inner solution is determined from Eq. (54) and reads

$\displaystyle
f_{b}^{0}(\tilde{r},\theta)=\underbrace{\frac{A_{0}}{\tilde{r}^{2}}\cos(\theta)}_{\text{must
match }F_{b}^{0}},f_{s}^{0}(\theta)=\tilde{f}_{s}^{0}\cos(\theta),\hskip
5.69046ptA_{0}=-\frac{\alpha}{2},\
\alpha=\frac{2D_{s}\xi(\gamma-1)+D_{b}}{D_{s}\xi(2\gamma+1)+D_{b}},\
\tilde{f}_{s}^{0}=\frac{\xi[2D_{s}\xi(2\gamma+1)-D_{b}]}{2[D_{s}\xi(2\gamma+1)+D_{b}]}.$
(55)

It us useful to introduce the pseudo-Green’s function $G$ with periodic
boundary conditions, which satisfies

$\displaystyle-\nabla^{2}G=\delta(\mathbf{r})-1,\ G(x\pm 1,y,z)=G(x,y\pm
1,z)=G(x,y,z\pm 1)=G(x,y,z).$ (56)

In fact, we will not need the full expression of $G$ but only its asymptotic
behavior near the origin, namely

$\displaystyle G(r,\theta,\phi)\underset{r\to 0}{=}\frac{1}{4\pi
r}+\frac{r^{2}}{6}+r^{4}A_{4}(\theta,\phi)+r^{8}A_{8}(\theta,\phi)+...$ (57)

with $A_{i}(\theta,\phi)$ functions that are built from Legendre polynomials
and respect the symmetry of the problem. For example, $A_{4}$ is built as a
combination $aP_{4}^{0}(\cos\theta)+bP_{4}^{4}(\cos\theta)\cos(4\phi)$ (with
$P_{n}^{m}$ the Legendre polynomials) with the requirement that, for
$\phi=n\pi/2$ (with $n$ an integer), $A_{4}$ does not depend on
$\cos(2\theta)$. The only possible form for $A_{4}$ is

$A_{4}(\theta,\phi)=C_{0}\left\\{20\cos(2\theta)+35\cos(4\theta)+9-5\cos(4\phi)[4\cos(2\theta)-\cos(4\theta)-3]\right\\},$
(58)

with $C_{0}$ a constant that we will not need to determine here. We claim that
$F_{b}^{0}(r,\theta)=-4\pi A_{0}\textbf{e}_{z}\cdot\nabla G$. Indeed, it is a
harmonic function with periodic boundary conditions, and its expansion near
the origin reads:

$F_{b}^{0}(r,\theta)=-4\pi A_{0}\textbf{e}_{z}\cdot\nabla
G=\underbrace{\frac{A_{0}}{r^{2}}\cos(\theta)}_{\text{ matching
}f_{b}^{0}}\underbrace{-\frac{4\pi}{3}A_{0}\cos(\theta)r}_{\text{must match
}f_{b}^{3}}+\underbrace{D_{0}[3\cos(\theta)+5\cos(3\theta)]r^{3}}_{\text{must
match }f_{b}^{5}}+\underbrace{E_{0}(\theta,\phi)r^{7}}_{\text{must match
}f_{b}^{9}}+o(r^{7}),$ (59)

with $D_{0}=-128\pi C_{0}$ and $E_{0}(\theta,\phi)$ a combination of Legendre
polynomials, for which we will only need the property:

$\int_{0}^{\pi}\mbox{d}\theta\int_{0}^{2\pi}\mbox{d}\phi\sin(\theta)\cos(\theta)E_{0}(\theta,\phi)=0.$
(60)

Next, at first and second order, we obtain the trivial solutions
$f_{b}^{1}=f_{s}^{1}=F_{b}^{1}=0$ and $f_{b}^{2}=f_{s}^{2}=F_{b}^{2}=0$. The
inner solution at third order is obtained by solving Eq. (54) with $i=3$ and
requiring that the large distance behavior matches with the corresponding term
in the expansion of $F_{b}^{0}$ in Eq. (59), we obtain

$\displaystyle
f_{b}^{3}(\tilde{r},\theta)=\underbrace{\frac{A_{3}}{\tilde{r}^{2}}\cos\theta}_{\text{must
match }F_{b}^{3}}+\underbrace{B_{3}\tilde{r}\cos\theta}_{\text{matching
}F_{b}^{0}},\ f_{s}^{3}(\theta)=\cos\theta\tilde{f}_{s}^{3},\
B_{3}=-\frac{4\pi}{3}A_{0},A_{3}=\frac{B_{3}}{2}\alpha,\
\tilde{f}_{s}^{3}=\frac{3\xi
D_{b}B_{3}}{2\left(D_{s}\xi(2\gamma+1)+D_{b}\right)}.$ (61)

The outer solution at this order and its expansion near the origin is

$F_{b}^{3}(r,\theta,\phi)=-\frac{4\pi}{3}A_{3}\textbf{e}_{z}\cdot\nabla
G=\underbrace{\frac{A_{3}}{r^{2}}\cos(\theta)}_{\text{ matching
}f_{b}^{3}}\underbrace{-\frac{4\pi}{3}A_{3}\cos(\theta)r}_{\text{must match
}f_{b}^{6}}+\underbrace{D_{1}[3\cos(\theta)+5\cos(3\theta)]r^{3}}_{\text{must
match }f_{b}^{8}}+\mathcal{O}(r^{7}).$ (62)

Next, the inner and outer solutions at fourth order vanish,
$f_{b}^{4}=f_{s}^{4}=F_{b}^{4}=0$. The inner solution $f_{b}^{5}$ reads:

$\displaystyle
f_{b}^{5}(\tilde{r},\theta)=\underbrace{\left[3\cos(\theta)+5\cos(3\theta)\right]\frac{A_{5}}{\tilde{r}^{4}}}_{\text{must
match
}F_{b}^{7}}+\underbrace{D_{0}\left[3\cos(\theta)+5\cos(3\theta)\right]\tilde{r}^{3}}_{\text{matching
}F_{b}^{0}}\ ,\
f_{s}^{5}(\theta)=\left[3\cos(\theta)+5\cos(3\theta)\right]\tilde{f}_{s}^{5}.$
(63)

The constants $A_{5}$ and $\tilde{f}_{s}^{5}$ are found by solving (54) with
$i=5$ but their precise expressions are not necessary because

$\int_{0}^{\pi}\mbox{d}\theta\int_{0}^{2\pi}\mbox{d}\phi\cos(\theta)\sin(\theta)\left[3\cos(\theta)+5\cos(3\theta)\right]=0,$
(64)

so that this order does not contribute in the expression of $D_{e}$. For the
outer solution, there are no terms to be matched so $F_{b}^{5}=0$. We continue
the iteration and we find at order $6$:

$\displaystyle
f_{b}^{6}=\underbrace{\frac{A_{6}}{\tilde{r}^{2}}\cos(\theta)}_{\text{must
match }F_{b}^{6}}+\underbrace{B_{6}\tilde{r}\cos(\theta)}_{\text{matching
}F_{b}^{3}}\ ,\ f_{s}^{6}(\theta)=\cos(\theta)\tilde{f}_{s}^{6},\
B_{6}=-\frac{4\pi}{3}A_{3},A_{6}=\frac{B_{6}}{2}\alpha,\
\tilde{f}_{s}^{6}=\frac{3\xi
D_{b}B_{6}}{2\left(D_{s}\xi(2\gamma+1)+D_{b}\right)},$ (65) $\displaystyle
F_{b}^{6}(r,\theta,\phi)=-4\pi A_{6}\textbf{e}_{z}\cdot\nabla
G=\underbrace{\frac{A_{6}}{r^{2}}\cos(\theta)}_{\text{ matching
}f_{b}^{6}}\underbrace{-\frac{4\pi}{3}A_{6}\cos(\theta)r}_{\text{must match
}f_{b}^{9}}+\underbrace{D_{2}[3\cos(\theta)+5\cos(3\theta)]r^{3}}_{\text{must
match }f_{b}^{11}}+\mathcal{O}(r^{7}).$ (66)

For the inner solution at $7^{th}$ order, there is no term to be matched, so
$f_{b}^{7}=f_{s}^{7}=0$. However, for the outer solution $F_{b}^{7}$, we have
to match a term from $f_{b}^{5}$. We can build $F_{b}^{7}$ from
$(\mathbf{e}_{z}\cdot\nabla)^{2}F_{b}^{0}$ so that its expansion near the
origin reads

$F_{b}^{7}(r,\theta)=\frac{4A_{5}}{3A_{0}}(\mathbf{e}_{z}\cdot\nabla)^{2}F_{b}^{0}=\underbrace{\frac{A_{5}}{r^{4}}\left[3\cos(\theta)+5\cos(3\theta)\right]}_{\text{
matching
}f_{b}^{5}}+\underbrace{\frac{64A_{5}D_{0}}{A_{0}}r\cos(\theta)}_{\text{must
match }f_{b}^{10}}+\mathcal{O}(r^{3}).$ (67)

$8^{th}$ order: Here, the form of the inner solution is

$\displaystyle
f_{b}^{8}(\tilde{r},\theta)=\underbrace{\left[3\cos(\theta)+5\cos(3\theta)\right]\frac{A_{8}}{\tilde{r}^{4}}}_{\text{must
match
}F_{b}^{10}}+\underbrace{\left[3\cos(\theta)+5\cos(3\theta)\right]D_{1}\tilde{r}^{3}}_{\text{matching
}F_{b}^{3}}\ ,\
f_{s}^{8}(\theta)=\tilde{f}_{s}^{8}\left[3\cos(\theta)+5\cos(3\theta)\right].$
(68)

As in the case of the $5^{th}$ order, these terms do not contribute to the
expression of $D_{e}$. There are no matching term for the outer solution so it
simply reads $F_{b}^{8}=0$. Finally, the inner solution at $9^{th}$ order
reads

$\displaystyle f_{b}^{9}(\tilde{r},\theta)$
$\displaystyle=\underbrace{\frac{A_{9}}{\tilde{r}^{2}}\cos(\theta)}_{\text{must
match }F_{b}^{9}}+\underbrace{B_{9}\tilde{r}\cos(\theta)}_{\text{matching
}F_{b}^{6}}+\underbrace{E_{0}(\theta,\phi)\tilde{r}^{7}}_{\text{matching
}F_{b}^{0}}+\underbrace{\frac{C_{9}E_{0}(\theta,\phi)}{\tilde{r}^{8}}}_{\text{must
match }F_{b}^{15}}\ ,\ $ (69) $\displaystyle f_{s}^{9}(\theta)$
$\displaystyle=\tilde{f}_{s}^{9}\cos(\theta),\ A_{9}=\frac{B_{9}}{2}\alpha,\
B_{9}=-\frac{4\pi}{3}A_{6},\ \tilde{f}_{s}^{9}=\frac{3\xi
D_{b}B_{9}}{2\left(D_{s}\xi(2\gamma+1)+D_{b}\right)}.$ (70)

Since $E_{0}$ satisfies Eq. (60), the two last terms in $f_{b}^{9}$ do not
contribute in the expression of $D_{e}$.

We stop here the iteration procedure and we consider the long-time dispersion
Eq. (51) which reads:

$\displaystyle
D_{e}=P_{b}^{st,0}\left\\{\\{D_{b}(1-\varphi)+2D_{s}\xi\varphi+D_{b}\varphi\left[A_{0}+\sum_{n=1}^{3}R^{3n}(A_{3n}+B_{3n})\right]-2D_{s}\varphi\sum_{n=0}^{3}R^{3n}\tilde{f}_{s}^{3n}\right\\}+o(R^{12}).$
(71)

Using the above results for $A_{n},B_{n},\tilde{f}_{s}^{n}$, this expression
can be simplified, leading to a re-summation and to the result

$D_{e}=D_{b}P_{b}^{st,0}\
\left\\{1-\frac{3\varphi\alpha}{2}\left[1-\frac{\alpha\varphi}{2}+\left(\frac{\alpha\varphi}{2}\right)^{2}-\left(\frac{\alpha\varphi}{2}\right)^{3}\right]\right\\}+o(\varphi^{4})=\frac{D_{b}}{1+(3\xi-1)\varphi}\
\frac{1-\alpha\varphi}{1+\frac{\alpha\varphi}{2}}+o(\varphi^{4}).$ (72)

###  D.3 Comparison with existing results

####  D.3.1 Monte Carlo data of Ref. [42]

Here, we compare the prediction of Eq. (15) to Monte Carlo simulations of Ref.
[42], as they appear in Figure 9 of their Supplemental Material, for periodic
arrays of weakly sticky crowders. The result of this comparison is shown in
Fig.2, and here we describe briefly how to identify the parameters of our
surface-mediated diffusion model with their short range interaction model. In
Ref. [42], the tracer particle is submitted to an exponentially decaying
potential $U(D)$, with $D$ the distance to the crowder’s wall. The link
between the parameter $\delta$ in our description and the potential $U(D)$ can
be identified by considering the simple one-dimensional situation in which the
tracer particle diffuses between $x=0$ and $x=H$ and is subjected to the
potential $U(x)$, assumed to vanish for distances larger than a length
$\lambda$, with $\lambda\ll H$. In this case, the probability $p_{\Sigma}$ to
find the particle at a distance from the surface less than $\varepsilon$ (with
$\lambda\ll\epsilon\ll H$) reads

$p_{\Sigma}=\frac{\int_{0}^{\epsilon}e^{-\beta
U(x)}\mbox{d}x}{\int_{0}^{H}e^{-\beta
U(x)}\mbox{d}x}=\frac{\int_{0}^{\epsilon}[e^{-\beta
U(x)}-1]\mbox{d}x+\varepsilon}{\int_{0}^{H}[e^{-\beta U(x)}-1]\mbox{d}x+H}\
\simeq\dfrac{\int_{0}^{\infty}\left(e^{-\beta
U(x)}-1\right)\mbox{d}x}{\int_{0}^{\infty}\left(e^{-\beta
U(x)}-1\right)\mbox{d}x+H}.$ (73)

In our formalism, we readily find $p_{\Sigma}=\delta/(\delta+H)$, so that the
expression of the length $\delta$ in terms of $U$ is

$\delta=\int_{0}^{\infty}\left(e^{-\beta U(D)}-1\right)\mbox{d}D.$ (74)

Next, since there is no potential barrier between the surface bound state and
the bulk, we are in the fast exchange limit, $k_{a},k_{d}\to\infty$, leading
to $\gamma=0$ in Eq. (15). Finally, Ref. [42] considers a tracer particle of
finite radius $r_{d}$ in presence of spherical crowders of radius $r_{c}$ and
volume fraction $\phi$, this is equivalent to considering a punctual tracer
particle in presence of crowders of radius $R=r_{c}+r_{d}$, and thus the
volume fraction in our description is $\varphi=\phi(1+r_{d}/r_{c})^{3}$. We
obtain $\varphi=0.57$. For this rather high value, the obstacles are actually
touching each other, but this is physical since the crowders are actually not
touching each other. It is actually very surprising that our formula can
correctly predicts $D_{e}$ for these parameters, as found in Fig. 4 (left) and
Fig. 2 in the main text.

####  D.3.2 Comparison with the result of macrotransport theory

Figure 4: Left: Comparison between Monte Carlo simulation data for dispersion
in 3D arrays of obstacles (circles: data of Fig. 9 in the Supplemental
Material of Ref. [42]), our theoretical prediction and the result (75) of Ref.
[51]. Parameters: $\varphi=0.57$, $D_{s}=D_{b}$, $\gamma=0$, and $\delta$ is
identified using Eq. (74) for an exponentially decaying potential $U(D)$.
Right. Comparison between the data obtained from the numerical resolution of
our formalism (squares), stochastic simulations (triangles) and analytical
predictions (see legend). Parameters: $\varphi=0.2$, $D_{s}=D_{b}$, and
$\gamma=0$ (blue) or $\gamma=0.1$ (orange).

In Ref. [51], using the formalism of macrotransport theory, a formula was
derived for the dispersion of tracers in a periodic array of spheres in a
model where tracers can diffuse inside the sphere, and in addition undergo
surface diffusion or diffusion in the bulk. The boundary conditions at the
interface in their approach is
$p_{s}(\mathbf{r}_{s},t)=k_{ed}p_{b}(\mathbf{r}_{s},t)$ and
$p_{\text{inside}}(\mathbf{r}_{s},t)=K_{ed}p_{b}(\mathbf{r}_{s},t)$ where
$k_{ed}$ and $K_{ed}$ are Henry’s laws constants (the suffix $ed$ refers to
notations of Ref. [51]), and $p_{\text{inside}}(\mathbf{r},t)$ is the
probability density of tracer position inside the spheres. In principle, if
one takes $K=0$ to prevent the tracers from entering into the obstacles, the
physical situation of Ref. [51] is identical to our formalism in the fast
exchange limit ($\gamma=0$). The result of Ref. [51] for a periodic array of
spheres in 3D is given by their equation (4.23) and reads with our notations:

$D_{e}=\frac{D_{b}}{1+(3\xi-1)\varphi}\left[1-\frac{3\varphi}{2\left(1+\frac{D_{s}}{D_{b}}\xi\right)+\varphi\left(1-2\frac{D_{s}}{D_{b}}\xi\right)}\right]+o(\varphi^{4})=\frac{D_{b}}{1+(3\xi-1)\varphi}\
\frac{1-\varphi}{1+\frac{\alpha\varphi}{2}}+o(\varphi^{4}),$ (75)

which is identified by writing $K_{ed}=\gamma_{ed}=0$,
$k_{ed}\alpha_{ed}=3\xi\varphi$, $\Gamma_{ed}=D_{s}\xi/D_{b}$. However, the
above equation (75) disagrees with our Eq. (15), since the numerator contains
a term $1-\varphi$ rather than $1-\alpha\varphi$.

In Fig. 4 (right), we show our prediction [Eq. (15)] along with the result
(75), compared with the numerical solution of our formalism. Our analytical
formula agrees with the numerical results, while the formula 75 does not and
in particular the non-monotonicity of $D_{e}$ with $\delta$ is absent. An
additional check of the validity of our formalism is provided by comparison
with direct stochastic numerical simulations of the process $\mathbf{r}(t)$
corresponding to the Fokker-Planck equations (1) and (2) (described in
Appendix C). The values of $D_{e}$ directly measured from the stochastic
trajectories obtained with this algorithm are shown on Fig. 4 (right) and
confirm the validity of our formalism.

##  Appendix E Dispersion in slowly undulated channels ($a\ll L$) - Proof of
Eq. (18)

###  E.1 The 2D symmetric channel

In this section, we determine the effective diffusivity $D_{e}$ in a two-
dimensional symmetric channel of height $h(x)$ satisfying $h(x+L)=h(x)$, in
the limit of large $L$. The channel wall is parametrized by $x$, we note that
$\mathbf{e}_{s,x}=\mathbf{e}_{x}+h^{\prime}(x)\mathbf{e}_{y}$, so that
$g_{xx}=1+h^{\prime}(x)^{2}$,
$\nabla_{s}^{2}x=-h^{\prime}(x)h^{\prime\prime}(x)/\left[1+h^{\prime}(x)^{2}\right]^{2}$.
Using the formulas of Appendix B, Eq. (12)-(13) read

$\left\\{\begin{split}&\frac{\partial^{2}f_{b}}{\partial
x^{2}}+\frac{\partial^{2}f_{b}}{\partial y^{2}}=0,\hskip
28.45274pt&0<y<h(x),\\\
&\frac{D_{s}}{\sqrt{1+h^{\prime}(x)^{2}}}\frac{\partial}{\partial
x}\left[\frac{1}{\sqrt{1+h^{\prime}(x)^{2}}}\frac{\partial f_{s}}{\partial
x}\right]-k_{d}f_{s}+k_{a}f_{b}=-D_{s}P_{s}^{st,0}\frac{h^{\prime}(x)h^{\prime\prime}(x)}{[1+h^{\prime}(x)^{2}]^{2}},\hskip
28.45274pt&y=h(x),\\\
&\frac{D_{b}}{\sqrt{1+h^{\prime}(x)^{2}}}\left[\frac{\partial f_{b}}{\partial
y}-h^{\prime}({x})\frac{\partial f_{b}}{\partial
x}\right]=k_{d}f_{s}-k_{a}f_{b}-\frac{D_{b}P_{b}^{st,0}h^{\prime}(x)}{\sqrt{1+h^{\prime}(x)^{2}}},\hskip
28.45274pt&y=h(x).\end{split}\right.$ (76)

Furthermore, $f_{b}$ and $f_{s}$ are periodic with period $L$, and satisfy
$\partial_{y}f_{b}|_{y=0}=0$ (from the symmetry $y\to-y$). The stationary
probability densities and the effective diffusivity are given by

$\displaystyle P_{b}^{st,0}=\frac{1}{V+S\delta},\hskip
28.45274ptP_{s}^{st,0}=\delta P_{b}^{st,0},\hskip 28.45274ptV=2\int_{0}^{L}dx\
h(x),\hskip 28.45274ptS=2\int_{0}^{L}dx\sqrt{1+h^{\prime}(x)^{2}},$ (77)
$\displaystyle
D_{e}=P_{b}^{st,0}VD_{b}+2D_{s}P_{s}^{st,0}\int_{0}^{L}\frac{dx}{\sqrt{1+h^{\prime}(x)^{2}}}-D_{b}\int_{0}^{L}dx\int_{-h(x)}^{h(x)}dy\frac{\partial
f_{b}}{\partial
x}-2D_{s}\int_{0}^{L}\frac{dx}{\sqrt{1+h^{\prime}(x)^{2}}}\frac{\partial
f_{s}}{\partial x}.$ (78)

We now analyze the above equations in the limit $L\to\infty$, when all other
parameters are held constant. Here we write $h(x)=\zeta(x/L)$, where the
function $\zeta(\tilde{x})$ does not depend on $L$, and we use the notation
$\tilde{x}=x/L$. The expansion of $f_{b}$ takes the form:

$\displaystyle f_{b}(\mathbf{r})=f_{b}^{0}(\tilde{x},y)+L^{-1}\
f_{b}^{1}(\tilde{x},y)+L^{-2}\ f_{b}^{2}(\tilde{x},y)+...,\hskip
14.22636ptf_{s}(x)=f_{s}^{0}(\tilde{x})+L^{-1}\ f_{s}^{1}(\tilde{x})+L^{-2}\
f_{s}^{2}(\tilde{x})+...$ (79)

Note that this expansion in powers of $L$ is equivalent to an expansion in
terms of the small parameter $\varepsilon=a/L$, with $a$ the minimal channel
width. This kind of calculation is very similar to those in lubrication theory
in hydrodynamics.

In order to perform the expansion it is useful to consider the equations in
terms of the rescaled variable $\tilde{x}=x/L$:

$\left\\{\begin{split}&\frac{1}{L^{2}}\frac{\partial^{2}f_{b}}{\partial\tilde{x}^{2}}+\frac{\partial^{2}f_{b}}{\partial
y^{2}}=0,&0<y<\zeta(\tilde{x}),\\\
&\frac{D_{s}}{1+\zeta^{\prime}(\tilde{x})^{2}/L^{2}}\left[\frac{1}{L^{2}}\frac{\partial^{2}\tilde{f}_{s}}{\partial\tilde{x}^{2}}-\frac{1}{L^{4}}\frac{\zeta^{\prime}(\tilde{x})\zeta^{\prime\prime}(\tilde{x})}{1+\zeta^{\prime}(\tilde{x})^{2}/L^{2}}\frac{\partial\tilde{f}_{s}}{\partial\tilde{x}}\right]=k_{d}f_{s}-k_{a}f_{b}-\frac{D_{s}P_{s}^{st,0}\zeta^{\prime}(\tilde{x})\zeta^{\prime\prime}(\tilde{x})}{L^{3}\left(1+\zeta^{\prime}(\tilde{x})^{2}/L^{2}\right)^{2}},&y=\zeta(\tilde{x})\\\
&\frac{D_{b}}{\sqrt{1+\zeta^{\prime}(\tilde{x})^{2}/L^{2}}}\left[\frac{\partial
f_{b}}{\partial y}-\frac{\zeta^{\prime}(\tilde{x})}{L^{2}}\frac{\partial
f_{b}}{\partial\tilde{x}}\right]\left[\tilde{x},h(\tilde{x})\right]=k_{d}f_{s}-k_{a}f_{b}-\frac{D_{b}P_{b}^{st,0}\zeta^{\prime}(\tilde{x})}{L\sqrt{1+\zeta^{\prime}(\tilde{x})^{2}/L^{2}}},&y=\zeta(\tilde{x}).\end{split}\right.$
(80)

with

$\displaystyle
P_{b}^{st,0}=\frac{1}{CL}\left(1+\mathcal{O}\left(L^{-2}\right)\right),\hskip
28.45274ptP_{s}^{st,0}=\frac{\delta}{CL}\left(1+\mathcal{O}\left(L^{-2}\right)\right),\hskip
28.45274ptC=2\left(\left<\zeta\right>+\delta\right).$ (81)

We insert the ansatz (79) into the above equations to find successive orders
$1/L^{n}$. At leading order, we obtain

$\displaystyle\frac{\partial^{2}f_{b}^{0}}{\partial y^{2}}=0,\hskip
28.45274pt(f_{s}^{0}-\delta f_{b}^{0})_{y=\zeta(\tilde{x})}=0,\hskip
28.45274pt\left(\frac{D_{b}}{k_{d}}\frac{\partial f_{b}^{0}}{\partial
y}-f_{s}^{0}+\delta f_{b}^{0}\right)_{y=\zeta(\tilde{x})}=0.$ (82)

Taking into account the condition $\partial_{y}f_{b}^{0}|_{y=0}=0$, the above
equations indicate that $f_{b}^{0}$ and $f_{s}^{0}$ do not depend on $y$, and
that they are proportional at the surface:

$f_{b}^{0}(\tilde{x},y)=f_{b}^{0}(\tilde{x})\ ,\ f_{s}^{0}(\tilde{x})=\delta
f_{b}^{0}(\tilde{x}).$ (83)

Next, at first order in $1/L$, we obtain the same equations as (83), leading
to $f_{b}^{1}(\tilde{x},y)=f_{b}^{1}(\tilde{x})$ and
$f_{s}^{1}(\tilde{x})=\delta f_{b}^{1}(\tilde{x}).$ At second order, we obtain

$\displaystyle\frac{\partial^{2}f_{b}^{0}}{\partial\tilde{x}^{2}}+\frac{\partial^{2}f_{b}^{2}}{\partial
y^{2}}$ $\displaystyle=0,$ (84) $\displaystyle
D_{s}\frac{\partial^{2}f_{s}^{0}}{\partial\tilde{x}^{2}}-k_{d}f_{s}^{2}+k_{a}f_{b}^{2}\
$ $\displaystyle,\ D_{b}\left[\frac{\partial f_{b}^{2}}{\partial
y}-\frac{\zeta^{\prime}(\tilde{x})^{2}}{2}\frac{\partial f_{b}^{0}}{\partial
y}-\zeta^{\prime}(\tilde{x})\frac{\partial
f_{b}^{0}}{\partial\tilde{x}}\right]=\left(k_{d}f_{s}^{2}-k_{a}f_{b}^{2}\right)-\frac{D_{b}}{C}\zeta^{\prime}(\tilde{x}),\hskip
28.45274pt[y=\zeta(\tilde{x})].$ (85)

Using the result at order $0$ [Eq. (83)], we see that Eq. (84) simplifies to

$\frac{\partial f_{b}^{2}}{\partial
y}(\tilde{x},y)=-{f_{b}^{0}}^{\prime\prime}(\tilde{x})y,$ (86)

where we have used the symmetry $y\to-y$. Next, combining Eqs. (85) and (86)
we obtain

${f_{b}^{0}}^{\prime\prime}(\tilde{x})\left(D_{b}\zeta(\tilde{x})+D_{s}\delta\right)+{f_{b}^{0}}^{\prime}(\tilde{x})D_{b}\zeta^{\prime}(\tilde{x})=\frac{D_{b}\zeta^{\prime}(\tilde{x})}{C}.$
(87)

The solution of this first-order differential equation that satisfies the
periodic boundary conditions is

${f_{b}^{0}}^{\prime}(\tilde{x})=\frac{1}{C\left(D_{b}\zeta(\tilde{x})+D_{s}\delta\right)}\left[D_{b}\zeta(\tilde{x})-\frac{\left<D_{b}\zeta/\left(D_{b}\zeta+D_{s}\delta\right)\right>}{\left<\left(D_{b}\zeta+D_{s}\delta\right)^{-1}\right>}\right].$
(88)

Inserting this result into Eq. (78), we find (after a few lines of algebra)

$D_{e}=\frac{1}{\left<h+\delta\right>\left<\left(D_{b}h+D_{s}\delta\right)^{-1}\right>}+\mathcal{O}(L^{-2}).$
(89)

###  E.2 The 3D axisymmetric channel

We now consider an axisymmetric channel oriented along the $x$ direction of
local radius $h(x)$. We denote by $r$ the distance to the central axis and
$\theta$ the angle around this axis. due to rotational invariance,
$f_{b}(r,x,\theta)$ and $f_{s}(x,\theta)$ do not depend on $\theta$. We
parametrize the surface with $x$ and $\theta$. The metric is identified by
calculating
$\textbf{e}_{s,x}=\partial_{x}\textbf{r}_{s}=\textbf{e}_{x}+h^{\prime}(x)\textbf{e}_{r}$
and $\textbf{e}_{s,\theta}=h(x)\textbf{e}_{\theta}$ with $\textbf{e}_{r}$ the
unit vector perpendicular to the channel axis, and
$\textbf{e}_{\theta}=\textbf{e}_{x}\times\textbf{e}_{r}$. We thus obtain
$g_{xx}=1+h^{\prime}(x)^{2}$, $g_{x,\theta}=0$ and
$g_{\theta\theta}=h(x)^{2}$. Using Appendix B, Eqs. (12,13) become

$\left\\{\begin{split}&\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial f_{b}}{\partial
r}\right)+\frac{\partial^{2}f_{b}}{\partial x^{2}}=0,\\\
&\frac{D_{s}}{1+h^{\prime}(x)^{2}}\left\\{\frac{\partial^{2}f_{s}}{\partial
x^{2}}+\left[\frac{h^{\prime}(x)}{h(x)}-\frac{h^{\prime}(x)h^{\prime\prime}(x)}{1+h^{\prime}(x)^{2}}\right]\frac{\partial
f_{s}}{\partial
x}\right\\}=k_{d}f_{s}-k_{a}f_{b}+\frac{D_{s}P_{s}^{st,0}}{1+h^{\prime}(x)^{2}}\left\\{\frac{h^{\prime}(x)}{h(x)}-\frac{h^{\prime}(x)h^{\prime\prime}(x)}{1+h^{\prime}(x)^{2}}\right\\},\\\
&\frac{D_{b}}{\sqrt{1+h^{\prime}(x)^{2}}}\left[\frac{\partial f_{b}}{\partial
r}-h^{\prime}(x)\frac{\partial f_{b}}{\partial
x}\right]\left[x,h(x)\right]=k_{d}f_{s}-k_{a}f_{b}-\frac{D_{b}P_{b}^{st,0}h^{\prime}(x)}{\sqrt{1+h^{\prime}(x)^{2}}},\end{split}\right.$
(90)

and $D_{e}$ is given by

$D_{e}=P_{b}^{st,0}VD_{b}+2\pi
D_{s}P_{s}^{st,0}\int_{0}^{L}\frac{dxh(x)}{\sqrt{1+h^{\prime}(x)^{2}}}-2\pi
D_{b}\int_{0}^{L}dx\int_{0}^{h(x)}dyy\frac{\partial f_{b}}{\partial x}-2\pi
D_{s}\int_{0}^{L}\frac{dxh(x)}{\sqrt{1+h^{\prime}(x)^{2}}}\frac{\partial
f_{s}}{\partial x}.$ (91)

As in the 2D case, we set $h(x)=\zeta(x/L)$ and we use the notation
$\tilde{x}=x/L$, in this notation we find

$\left\\{\begin{split}&\frac{1}{r}\frac{\partial}{\partial
r}\left(r\frac{\partial f_{b}}{\partial
r}\right)+\frac{1}{L^{2}}\frac{\partial^{2}f_{b}}{\partial\tilde{x}^{2}}=0,\\\
&\frac{D_{s}}{L^{2}+\zeta^{\prime}(\tilde{x})^{2}}\left\\{\frac{\partial^{2}f_{s}}{\partial\tilde{x}^{2}}+\left[\frac{\zeta^{\prime}(\tilde{x})}{\zeta(\tilde{x})}-\frac{\zeta^{\prime}(\tilde{x})\zeta^{\prime\prime}(\tilde{x})}{L^{2}+\zeta^{\prime}(\tilde{x})^{2}}\right]\frac{\partial
f_{s}}{\partial\tilde{x}}\right\\}=k_{d}f_{s}-k_{a}f_{b}+\frac{LD_{s}P_{s}^{st,0}}{L^{2}+\zeta^{\prime}(\tilde{x})^{2}}\left\\{\frac{\zeta^{\prime}(\tilde{x})}{\zeta(\tilde{x})}-\frac{\zeta^{\prime}(\tilde{x})\zeta^{\prime\prime}(\tilde{x})}{\left[L^{2}+\zeta^{\prime}(\tilde{x})^{2}\right]}\right\\},\\\
&\frac{D_{b}}{\sqrt{1+\zeta^{\prime}(\tilde{x})^{2}/L^{2}}}\left[\frac{\partial
f_{b}}{\partial r}-\frac{\zeta^{\prime}(\tilde{x})}{L^{2}}\frac{\partial
f_{b}}{\partial\tilde{x}}\right]_{\left[\tilde{x},\zeta(\tilde{x})\right]}=k_{d}f_{s}-k_{a}f_{b}-\frac{D_{b}P_{b}^{st,0}\zeta^{\prime}(\tilde{x})}{L\sqrt{1+\zeta^{\prime}(\tilde{x})^{2}/L^{2}}},\end{split}\right.$
(92)

with

$\displaystyle
P_{b}^{st,0}=\frac{1}{C^{\prime}L}\left[1+\mathcal{O}\left(L^{-2}\right)\right],\
P_{s}^{st,0}=\frac{\delta}{C^{\prime}L}\left[1+\mathcal{O}\left(L^{-2}\right)\right],\
C^{\prime}=\pi\left<h^{2}\right>+2\pi\delta\left<h\right>.$ (93)

In the limit $L\to\infty$, the solutions have a series expansion (see Eq.
(17)) and successive orders are found by applying the same method as for the
2D channel. As in the 2D case, at leading and first order, we find

$f_{b}^{0}(\tilde{x},r)=f_{b}^{0}(\tilde{x}),\ f_{s}^{0}(\tilde{x})=\delta
f_{b}^{0}(\tilde{x}),\ f_{b}^{1}(\tilde{x},r)=f_{b}^{1}(\tilde{x}),\
f_{s}^{1}(\tilde{x})=\delta f_{b}^{1}(\tilde{x}).$ (94)

The equations at $2^{nd}$ order are

$\displaystyle\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial
f_{b}^{2}}{\partial
r}\right)+\frac{\partial^{2}f_{b}^{0}}{\partial\tilde{x}^{2}}=0,$ (95)
$\displaystyle
D_{s}\left[\frac{\partial^{2}f_{s}^{0}}{\partial\tilde{x}^{2}}+\frac{\zeta^{\prime}(\tilde{x})}{\zeta(\tilde{x})}\frac{\partial
f_{s}^{0}}{\partial\tilde{x}}\right]_{\left[\tilde{x},\zeta(\tilde{x})\right]}=k_{d}f_{s}^{2}-k_{a}f_{b}^{2}+\frac{D_{s}\delta\zeta^{\prime}(\tilde{x})}{C^{\prime}\zeta(\tilde{x})},$
(96) $\displaystyle D_{b}\left[\frac{\partial f_{b}^{2}}{\partial
r}-\frac{\zeta^{\prime}(\tilde{x})^{2}}{2}\frac{\partial f_{b}^{0}}{\partial
r}-\zeta^{\prime}(\tilde{x})\frac{\partial
f_{b}^{0}}{\partial\tilde{x}}\right]_{\left[\tilde{x},\zeta(\tilde{x})\right]}=k_{d}f_{s}^{2}-k_{a}f_{b}^{2}-\frac{D_{b}\zeta^{\prime}(\tilde{x})}{C^{\prime}}.$
(97)

The integration of Eq. (95) leads to
$\partial_{r}f_{b}^{2}=-\frac{r{f_{b}^{0}}^{\prime\prime}(\tilde{x})}{2}$
which can be combined with Eqs. (96, 97), leading to

${f_{b}^{0}}^{\prime\prime}(\tilde{x})\left[\frac{D_{b}\zeta^{2}(\tilde{x})}{2}+D_{s}\delta\zeta(\tilde{x})\right]+{f_{b}^{0}}^{\prime}(\tilde{x})\left[D_{b}\zeta(x)\zeta^{\prime}(\tilde{x})+D_{s}\delta\zeta^{\prime}(\tilde{x})\right]=\frac{1}{C^{\prime}}\left[D_{b}\zeta(\tilde{x})\zeta^{\prime}(\tilde{x})+D_{s}\delta\zeta^{\prime}(\tilde{x})\right].$
(98)

Solving this differential equation, we obtain:

${f_{b}^{0}}^{\prime}(\tilde{x})=\frac{1}{C^{\prime}}\left[1-\frac{1}{\zeta(\tilde{x})\left(D_{b}\zeta(\tilde{x})+2D_{s}\delta\right)\left<h^{-1}\left(D_{b}h+2D_{s}\delta\right)^{-1}\right>}\right].$
(99)

Inserting this result into Eq. (91) yields the $D_{e}$ at leading order in
$1/L$:

$D_{e}=\frac{1}{\left<h\left(h+2\delta\right)\right>\left<h^{-1}\left(D_{b}h+2D_{s}\delta\right)^{-1}\right>}+\mathcal{O}\left(L^{-2}\right).$
(100)

###  E.3 Upper bounds on critical surface diffusion $D_{s}^{*}$

Here we show that

$D_{s}^{*}\leq D_{s}^{0}\leq D_{b},$ (101)

where we recall that $D_{s}^{*}$ is the value of $D_{s}$ for which $\partial
D_{e}/\partial\delta=0$, and $D_{s}^{0}$ is the value of $D_{s}$ for which
$D_{e}(\delta=0)=D_{e}(\delta=\infty)$. Using Eq. (18) it is straight forward
to show that

$\displaystyle
D_{s}^{*}=D_{b}\frac{\left<h^{d-2}\right>\left<h^{1-d}\right>}{\left<h^{d-1}\right>\left<h^{-d}\right>},$
$\displaystyle\
D_{s}^{0}=D_{b}\frac{\left<h^{d-2}\right>\left<h^{2-d}\right>}{\left<h^{d-1}\right>\left<h^{1-d}\right>}.$
(102)

Here we apply the Cauchy-Schwarz inequality which states that, for all
integrable functions $f,g$:

$\displaystyle\left<fg\right>^{2}\leq\left<f^{2}\right>\left<g^{2}\right>,$
(103)

The choice $f=h^{1-d/2}$ and $g=h^{-d/2}$ leads to the inequality
$\left<h^{1-d}\right>^{2}\leq\left<h^{-d}\right>\left<h^{2-d}\right>,$ which
leads to $D_{s}^{*}\leq D_{s}^{0}$. Now, using again Eq. (103), we find
$\left<h^{d-2}\right>^{2}\leq\left<h^{d-1}\right>\left<h^{d-3}\right>$ and
$\left<h^{2-d}\right>^{2}\leq\left<h^{1-d}\right>\left<h^{3-d}\right>$,
leading to

$\displaystyle\frac{\left<h^{d-2}\right>\left<h^{2-d}\right>}{\left<h^{d-1}\right>\left<h^{1-d}\right>}\leq\frac{\left<h^{d-3}\right>\left<h^{3-d}\right>}{\left<h^{d-2}\right>\left<h^{2-d}\right>}.$
(104)

By induction, we find

$\displaystyle\frac{D_{s}^{0}}{D_{b}}=\frac{\left<h^{d-2}\right>\left<h^{2-d}\right>}{\left<h^{d-1}\right>\left<h^{1-d}\right>}\leq\frac{\left<h^{d-3}\right>\left<h^{3-d}\right>}{\left<h^{d-2}\right>\left<h^{2-d}\right>}\leq...\leq\frac{1}{\left<h\right>\left<h^{-1}\right>}\leq
1,$ (105)

using $\left<h\right>\left<h^{-1}\right>\geq 1$. Thus, we recover Eq. (101).

###  E.4 Physical derivation for dispersion in slowly varying channels

Here, we see how the above results can be obtained via a physical argument
similar to that used for the original derivation of the Fick-Jacobs
approximation. For convenience, we consider a channel with a finite size
surface region of thickness $\ell$ where there is a uniform potential $v_{s}$,
in terms of the coordinates perpendicular to the channel axis the overall
potential $V({\bf r}_{\perp})$ is taken to be zero outside the surface region.
As the surface varies slowly along the channel, we assume that the position of
the particle in the height direction (the plane in the channel at longitudinal
position $x$ denoted $C(x)$) equilibrates giving an $x$ dependent partition
function

$Z(x)=\int_{C(x)}d{\bf r}_{\perp}e^{-\beta V({\bf
r}_{\perp})}=s_{d-1}h(x)^{d-1}+a_{d-1}h(x)^{d-2}\ell\exp(-\beta v_{s}),$ (106)

where $s_{d}$ denotes the volume of a unit $d$ dimensional sphere and $a_{d}$
its area. To recover the surface mediated diffusion model considered here one
takes the limit $\ell\to 0$ with $\ell\exp(-\beta v_{s})=\delta$ [see Eq.
(74)], and we note that $a_{d}=ds_{d}$. The free energy $\mathcal{F}$ as a
function of $x$ is thus given by

$\mathcal{F}(x)=-\frac{1}{\beta}\left[\ln\left(h(x)^{d-1}+(d-1)\delta
h(x)^{d-2}\right)+\ln(s_{d-1})\right].$ (107)

Also, we can introduce the average diffusion constant at position $x$ given by

$D_{a}(x)=\frac{D_{s}(d-1)\delta h(x)^{d-2}+D_{b}h(x)^{d-1}}{(d-1)\delta
h(x)^{d-2}+h(x)^{d-1}}.$ (108)

The effective diffusion constant in one dimension in a periodic potential
$\mathcal{F}(x)$ with periodic diffusion constant $D_{a}(x)$ is given by the
Lifson-Jackson formula [52]

$D_{e}=\frac{1}{\langle\exp(-\beta\mathcal{F})\rangle\langle\exp(\beta\mathcal{F})D^{-1}_{a}\rangle}=\frac{1}{\langle
h^{d-1}+(d-1)\delta h^{d-2}\rangle\langle\left(D_{b}h^{d-1}+D_{s}(d-1)\delta
h^{d-2}\right)^{-1}\rangle},$ (109)

which is precisely Eq. (18) of our Letter.

##  Appendix F Dispersion in highly undulated channels ($a\gg L$) - Proof of
Eq. (20)

This last section concerns channels in the large corrugation limit. Let us
call $r$ the distance to the central axis, and $w(r)$ the “half-width” of the
corrugations of the channel, see Fig. 5. For simplicity, we assume that the
corrugation is symmetric about x = 0. This assumption can in fact be be
avoided by slightly adapting the arguments presented below and does not change
the final result. In 2D, the surface is parametrized with the variable $r$.
The metric is determined by calculating
$\mathbf{e}_{s,r}=\mathbf{e}_{r}+w^{\prime}(r)\mathbf{e}_{x}$ which leads to
$g_{rr}=1+w^{\prime 2}(r)$. For an axisymmetric 3D channel, we parametrize the
surface with $r$ and $\theta$, the components of the metric read
$g_{rr}=1+w^{\prime 2}(r)$, $g_{r\theta}=0$ and $g_{\theta\theta}=r^{2}$.
Here, we consider the large corrugation limit $a\to\infty$, by assuming that
$w(r)=\tilde{w}(r/a)$, where $\tilde{w}$ does not depend on $a$. In terms of
the rescaled variable $\tilde{r}=r/a$, the equations read (for $d=2$ or
$d=3$):

$\left\\{\begin{split}&\text{(bulk) }\frac{\partial^{2}\tilde{f}_{b}}{\partial
x^{2}}+\frac{1}{a^{2}\tilde{r}^{d-2}}\frac{\partial}{\partial\tilde{r}}\left(\tilde{r}^{d-2}\frac{\partial\tilde{f}_{b}}{\partial\tilde{r}}\right)=0,\\\
&\text{(surface) }\frac{D_{s}}{a^{2}+\tilde{w}^{\prime
2}(\tilde{r})}\left\\{\frac{\partial^{2}\tilde{f}_{s}}{\partial\tilde{r}^{2}}-\delta\tilde{w}^{\prime\prime}(\tilde{r})+\left[\frac{d-2}{\tilde{r}}-\frac{\tilde{w}^{\prime}(\tilde{r})\tilde{w}^{\prime\prime}(\tilde{r})}{a^{2}+\tilde{w}^{\prime
2}(\tilde{r})}\right]\left[\frac{\partial\tilde{f}_{s}}{\partial\tilde{r}}-\delta\tilde{w}^{\prime}(\tilde{r})\right]\right\\}=k_{d}\tilde{f}_{s}-k_{a}\tilde{f}_{b},\\\
&\text{(b. c.) }\frac{D_{b}}{\sqrt{1+\frac{\tilde{w}^{\prime
2}(\tilde{r})}{a^{2}}}}\left[-\frac{\tilde{w}^{\prime}(\tilde{r})}{a^{2}}\frac{\partial\tilde{f}_{b}}{\partial\tilde{r}}+\frac{\partial\tilde{f}_{b}}{\partial
x}\right]\left[\tilde{w}(\tilde{r}),\tilde{r}\right]=k_{d}\tilde{f}_{s}-k_{a}\tilde{f}_{b}+\frac{D_{b}}{\sqrt{1+\frac{\tilde{w}^{\prime
2}(\tilde{r})}{a^{2}}}},\end{split}\right.$ (110)

where we used the notation $\tilde{f}_{b}=f_{b}/P_{b,0}^{st}$ and
$\tilde{f}_{s}=f_{s}/P_{b,0}^{st}$. At leading order in the limit
$a\to\infty$, we find

$\begin{split}\frac{\partial^{2}\tilde{f}_{b}^{0}}{\partial x^{2}}=0,\hskip
28.45274pt\left(k_{d}\tilde{f}_{s}^{0}-k_{a}\tilde{f}_{b}^{0}\right)_{\left[\tilde{w}(\tilde{r}),\tilde{r}\right]}=0,\hskip
28.45274ptD_{b}\left(\frac{\partial\tilde{f}_{b}^{0}}{\partial
x}\right)_{\left[\tilde{w}(\tilde{r}),\tilde{r}\right]}=k_{a}\tilde{f}_{b}^{0}-k_{d}\tilde{f}_{s}^{0}+D_{b}.\end{split}$
(111)

Solving these equations leads to:

$f_{b}^{0}(x,\tilde{r})=\begin{cases}A^{C}(\tilde{r})&\text{ if
}|\tilde{r}|<1,\\\ A^{P}(\tilde{r})+x&\text{ if }|\tilde{r}|\geq
1,\end{cases}$ (112)

where the functions $A^{C}$ and $A^{P}$ depend only on $\tilde{r}$. The first
line of the above expression is the solution in the central region which
respects periodic boundary conditions. The second one corresponds to the
peripheral region and satisfies the boundary condition given in (111). The
latter gives a contribution in the expression of the long-time dispersion at
leading order in $1/a$. Inserting the solution for $f$ into

$D_{e}\underset{a\rightarrow+\infty}{=}D_{b}P_{b,0}^{st}\left(V-V_{p}\right)=D_{b}P_{b,0}^{st}V_{c},$
(113)

with $V_{p}$ and $V_{c}$ the volumes of the peripheral and central region
respectively. Using $P_{b,0}^{st}=\left[V+\delta S\right]^{-1}$, we thus
recover Eq. (20).

Figure 5: Schematic drawing of a channel where the half-width $w$ is shown.

## References

  * [1] Marbach, S., Dean, D. S. & Bocquet, L.  Transport and dispersion across wiggling nanopores.  _Nat. Phys._ 14, 1108 (2018).
  * [2] Aminian, M., Bernardi, F., Camassa, R., Harris, D. M. & McLaughlin, R. M.  How boundaries shape chemical delivery in microfluidics.  _Science_ 0532 (2016).
  * [3] Kim, W. K., Kanduč, M., Roa, R. & Dzubiella, J.  Tuning the permeability of dense membranes by shaping nanoscale potentials.  _Phys. Rev. Lett._ 122, 108001 (2019).
  * [4] Brenner, H. & Edwards, D. A.  _Macrotransport processes_ (Butterworth-Heinemann, Boston, 1993).
  * [5] Le Borgne, T., Dentz, M. & Villermaux, E.  Stretching, coalescence, and mixing in porous media.  _Phys. Rev. Lett._ 110, 204501 (2013).
  * [6] Dentz, M., Le Borgne, T., Englert, A. & Bijeljic, B.  Mixing, spreading and reaction in heterogeneous media: A brief review.  _J. Contam. Hydrol._ 120, 1–17 (2011).
  * [7] Barros, F. P., Dentz, M., Koch, J. & Nowak, W.  Flow topology and scalar mixing in spatially heterogeneous flow fields.  _Geophys. Res. Lett._ 39 (2012).
  * [8] Bernate, J. A. & Drazer, G.  Stochastic and deterministic vector chromatography of suspended particles in one-dimensional periodic potentials.  _Phys. Rev. Lett._ 108, 214501 (2012).
  * [9] Aminian, M., Bernardi, F., Camassa, R. & McLaughlin, R. M.  Squaring the circle: Geometric skewness and symmetry breaking for passive scalar transport in ducts and pipes.  _Phys. Rev. Lett._ 115, 154503 (2015).
  * [10] Taylor, G.  Dispersion of soluble matter in solvent flowing slowly through a tube.  _Proc. R. Soc. Lon. A_ 219, 186–203 (1953).
  * [11] Dean, D. S., Drummond, I. & Horgan, R.  Effective transport properties for diffusion in random media.  _J. Stat. Mech.: Theor. Exp._ 2007, P07013 (2007).
  * [12] Bhattacharjee, T. & Datta, S. S.  Bacterial hopping and trapping in porous media.  _Nature communications_ 10, 1–9 (2019).
  * [13] Carusela, M. F. & Miguel Rubi, J.  Computational model for membrane transporters. potential implications for cancer.  _Frontiers in Cell and Developmental Biology_ 9, 333 (2021).
  * [14] Malgaretti, P., Pagonabarraga, I. & Rubi, M.  Entropic transport in confined media: a challenge for computational studies in biological and soft-matter systems.  _Frontiers in Physics_ 1, 21 (2013).
  * [15] Burada, P. S., Hänggi, P., Marchesoni, F., Schmid, G. & Talkner, P.  Diffusion in confined geometries.  _ChemPhysChem_ 10, 45–54 (2009).
  * [16] Jacobs, M.  _Diffusion processes_ (Springer, New-York, 1967).
  * [17] Reguera, D., Schmid, G., Burada, P. S., Rubi, J. M., Reimann, P., Hänggi, P.,  Entropic transport: Kinetics, scaling, and control mechanisms.  _Phys. Rev. Lett._ 96, 130603 (2006).
  * [18] Rubi, J. M.  Entropic diffusion in confined soft-matter and biological systems.  _Europhys. Lett._ 127, 10001 (2019).
  * [19] Yang, X. Liu, C. Li, Y. Marchesoni, F. Hänggi, P. Zhang, H. P.  Hydrodynamic and entropic effects on colloidal diffusion in corrugated channels.  _Proc. Natl. Acad. Sci. U. S. A._ 201707815 (2017).
  * [20] Martens, S., Straube, A., Schmid, G., Schimansky-Geier, L. & Hänggi, P.  Hydrodynamically enforced entropic trapping of brownian particles.  _Phys. Rev. Lett._ 110, 010601 (2013).
  * [21] Kalinay, P. & Percus, J.  Corrections to the fick-jacobs equation.  _Phys. Rev. E_ 74, 041203 (2006).
  * [22] Kalinay, P.  Nonscaling calculation of the effective diffusion coefficient in periodic channels.  _J. Chem. Phys._ 146, 034109 (2017).
  * [23] Martens, S., Schmid, G., Schimansky-Geier, L. & Hänggi, P.  Entropic particle transport: Higher-order corrections to the fick-jacobs diffusion equation.  _Phys. Rev. E_ 83, 051135 (2011).
  * [24] Zwanzig, R.  Diffusion past an entropy barrier.  _J Phys. Chem._ 96, 3926–3930 (1992).
  * [25] Mangeat, M., Guerin, T. & Dean, D. S.  Dispersion in two dimensional channels—the fick–jacobs approximation revisited.  _J. Stat. Mech.: Theor. Exp._ 2017, 123205 (2017).
  * [26] Mangeat, M., Guérin, T. & Dean, D. S.  Geometry controlled dispersion in periodic corrugated channels.  _Europhys. Lett._ 118, 40004 (2017).
  * [27] Mangeat, M., Guérin, T. & Dean, D.  Dispersion in two-dimensional periodic channels with discontinuous profiles.  _J. Chem. Phys._ 149, 124105 (2018).
  * [28] Brenner, H.  The slow motion of a sphere through a viscous fluid towards a plane surface.  _Chem. Eng. Sci._ 16, 242–251 (1961).
  * [29] Bychuk, O. V. & O’Shaughnessy, B.  Anomalous diffusion at liquid surfaces.  _Physical review letters_ 74, 1795 (1995).
  * [30] Walder, R., Nelson, N. & Schwartz, D. K.  Single molecule observations of desorption-mediated diffusion at the solid-liquid interface.  _Phys. Rev. Lett._ 107, 156102 (2011).
  * [31] Skaug, M. J., Mabry, J. & Schwartz, D. K.  Intermittent molecular hopping at the solid-liquid interface.  _Phys. Rev. Lett._ 110, 256101 (2013).
  * [32] Chee, S. W., Baraissov, Z., Loh, N. D., Matsudaira, P. T. & Mirsaidov, U.  Desorption-mediated motion of nanoparticles at the liquid–solid interface.  _J. Phys. Chem. C_ 120, 20462–20470 (2016).
  * [33] Wang, D. & Schwartz, D. K.  Non-brownian interfacial diffusion: Flying, hopping, and crawling.  _J. Phys. Chem. C_ 124, 19880–19891 (2020).
  * [34] Morrin, G. T., Kienle, D. F., Weltz, J. S., Traeger, J. C. & Schwartz, D. K.  Polyelectrolyte surface diffusion in a nanoslit geometry.  _Macromolecules_ 53, 4110–4120 (2020).
  * [35] Bénichou, O., Grebenkov, D., Levitz, P., Loverdo, C. & Voituriez, R.  Optimal reaction time for surface-mediated diffusion.  _Phys. Rev. Lett._ 105, 150606 (2010).
  * [36] Calandre, T., Bénichou, O. & Voituriez, R.  Accelerating search kinetics by following boundaries.  _Phys. Rev. Lett._ 112, 230601 (2014).
  * [37] Rupprecht, J.-F., Bénichou, O., Grebenkov, D. & Voituriez, R.  Exact mean exit time for surface-mediated diffusion.  _Phys. Rev. E_ 86, 041135 (2012).
  * [38] Monserud, J. H. & Schwartz, D. K.  Interfacial molecular searching using forager dynamics.  _Phys. Rev. Lett._ 116, 098303 (2016).
  * [39] von Hippel, P. H. & Berg, O. G.  Facilitated target location in biological systems.  _J. Biol. Chem._ 264, 675–678 (1989).
  * [40] Coppey, M., Bénichou, O., Voituriez, R. & Moreau, M.  Kinetics of target site localization of a protein on dna: a stochastic approach.  _Biophys. J._ 87, 1640–1649 (2004).
  * [41] Berg, O. G., Winter, R. B. & Von Hippel, P. H.  Diffusion-driven mechanisms of protein translocation on nucleic acids. 1. models and theory.  _Biochemistry_ 20, 6929–6948 (1981).
  * [42] Putzel, G. G., Tagliazucchi, M. & Szleifer, I.  Nonmonotonic diffusion of particles among larger attractive crowding spheres.  _Phys. Rev. Lett._ 113, 138302 (2014).
  * [43] Levesque, M., Bénichou, O., Voituriez, R. & Rotenberg, B.  Taylor dispersion with adsorption and desorption.  _Phys. Rev. E_ 86, 036316 (2012).
  * [44] Berezhkovskii, A. M. & Skvortsov, A. T.  Aris-taylor dispersion with drift and diffusion of particles on the tube wall.  _J. Chem. Phys._ 139, 084101 (2013).
  * [45] Quintard, M. & Whitaker, S.  Convection, dispersion, and interfacial transport of contaminants: Homogeneous porous media.  _Adv. Water. Ress._ 17, 221–239 (1994).
  * [46] Edwards, D. A., Shapiro, M. & Brenner, H.  Dispersion and reaction in two-dimensional model porous media.  _Physics of Fluids A: Fluid Dynamics_ 5, 837–848 (1993).
  * [47] Grebenkov, D. S.  Imperfect diffusion-controlled reactions.  _Chemical Kinetics: Beyond the Textbook_ 191–219 (2019).
  * [48] Alexandre, A., Guérin, T. & Dean, D. S.  Generalized taylor dispersion for translationally invariant microfluidic systems.  _Phys. Fluids_ 33, 082004 (2021).
  * [49] Guérin, T. & Dean, D. S.  Kubo formulas for dispersion in heterogeneous periodic nonequilibrium systems.  _Phys. Rev. E_ 92, 062103 (2015).
  * [50] Guérin, T. & Dean, D. S.  Force-induced dispersion in heterogeneous media.  _Phys. Rev. Lett._ 115, 020601 (2015).
  * [51] Edwards, D. A. & Davis, A.  Diffusion and convective dispersion through arrays of spheres with surface adsorption, diffusion, and unequal solute partitioning.  _Chemical engineering science_ 50, 1441–1454 (1995).
  * [52] Lifson, S. & Jackson, J. L.  On the Self-Diffusion of Ions in a Polyelectrolyte Solution.  _J. Chem. Phys._ 36, 2410–2414 (1962).
  * [53] Singer, A. & Schuss, Z. & Osipov, A. & Holcman, D.  Partially reflected diffusion.  _SIAM J. Appl. Math._ 68, 844–868 (2008).
  * [54] Ward, M. J. & Keller, J. B.  Strong localized perturbations of eigenvalue problems.  _SIAM J. Appl. Math._ 53, 770-798 (1993).
  * [55] Mamode, M.  Fundamental solution of the Laplacian on flat tori and boundary value problems for the planar Poisson equation in rectangles.  _Boundary Value Problems_ , 1-9 (2014).