A Alexandre

Thomas Guérin

M Mangeat

D S Dean

How stickiness can speed up diffusion in confined systems

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 [START_REF] Marbach | Transport and dispersion across wiggling nanopores[END_REF][START_REF] Aminian | How boundaries shape chemical delivery in microfluidics[END_REF][START_REF] Kim | Tuning the permeability of dense membranes by shaping nanoscale potentials[END_REF][START_REF] Brenner | Macrotransport processes[END_REF]. The effective diffusivity is a crucial input for problems of mixing [START_REF] Le Borgne | Stretching, coalescence, and mixing in porous media[END_REF][START_REF] Dentz | Mixing, spreading and reaction in heterogeneous media: A brief review[END_REF][START_REF] Barros | Flow topology and scalar mixing in spatially heterogeneous flow fields[END_REF], sorting [START_REF] Bernate | Stochastic and deterministic vector chromatography of suspended particles in one-dimensional periodic potentials[END_REF], chemical delivery [START_REF] Aminian | How boundaries shape chemical delivery in microfluidics[END_REF][START_REF] Aminian | Squaring the circle: Geometric skewness and symmetry breaking for passive scalar transport in ducts and pipes[END_REF] as well as chemical reactions [START_REF] Brenner | Macrotransport processes[END_REF][START_REF] Dentz | Mixing, spreading and reaction in heterogeneous media: A brief review[END_REF]. 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 [START_REF] Taylor | Dispersion of soluble matter in solvent flowing slowly through a tube[END_REF], and the decrease of dispersion due to the energy barriers created by time independent potentials [START_REF] Dean | Effective transport properties for diffusion in random media[END_REF]. In a number of systems, such as porous media [START_REF] Bhattacharjee | Bacterial hopping and trapping in porous media[END_REF], zeolites, and biological channels [START_REF] Carusela | Computational model for membrane transporters. potential implications for cancer[END_REF], 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 [START_REF] Malgaretti | Entropic transport in confined media: a challenge for computational studies in biological and soft-matter systems[END_REF][START_REF] Burada | Diffusion in confined geometries[END_REF][START_REF] Jacobs | Diffusion processes[END_REF][START_REF] Reguera | Entropic transport: Kinetics, scaling, and control mechanisms[END_REF][START_REF] Rubi | Entropic diffusion in confined soft-matter and biological systems[END_REF]. The effect of the confining geometry on effective diffusivity has been widely studied [START_REF] Aminian | How boundaries shape chemical delivery in microfluidics[END_REF][START_REF] Jacobs | Diffusion processes[END_REF][START_REF] Yang | Hydrodynamic and entropic effects on colloidal diffusion in corrugated channels[END_REF][START_REF] Martens | Hydrodynamically enforced entropic trapping of brownian particles[END_REF][START_REF] Kalinay | Corrections to the fick-jacobs equation[END_REF][START_REF] Kalinay | Nonscaling calculation of the effective diffusion coefficient in periodic channels[END_REF][START_REF] Martens | Entropic particle transport: Higher-order corrections to the fick-jacobs diffusion equation[END_REF][START_REF] Zwanzig | Diffusion past an entropy barrier[END_REF][START_REF] Mangeat | Dispersion in two dimensional channels-the Fick-Jacobs approximation revisited[END_REF][START_REF] Mangeat | Geometry controlled dispersion in periodic corrugated channels[END_REF][START_REF] Mangeat | Dispersion in two-dimensional periodic channels with discontinuous profiles[END_REF], 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 than that of simple Brownian motion at a reflecting wall. First, the tracer will typically be subject to non-specific interactions with the surface, for example as the result of Van der Waals long range attraction and electrostatic potentials, that can be attractive or repulsive, depending on geometry and charges [START_REF] Kim | Tuning the permeability of dense membranes by shaping nanoscale potentials[END_REF][START_REF] Israelachvili | Intermolecular and surface forces[END_REF]. 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 [START_REF] Brenner | The slow motion of a sphere through a viscous fluid towards a plane surface[END_REF]. 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 [START_REF] Bychuk | Anomalous diffusion at liquid surfaces[END_REF]. The stochastic motion resulting from interplay of bulk and surface diffusion has recently been directly observed experimentally [START_REF] Walder | Single molecule observations of desorption-mediated diffusion at the solid-liquid interface[END_REF][START_REF] Skaug | Intermittent molecular hopping at the solid-liquid interface[END_REF][START_REF] Chee | Desorption-mediated motion of nanoparticles at the liquid-solid interface[END_REF][START_REF] Wang | Non-brownian interfacial diffusion: Flying, hopping, and crawling[END_REF][START_REF] Morrin | Polyelectrolyte surface diffusion in a nanoslit geometry[END_REF].

It is also known that surface mediated transport gives rise to non-trivial effects in the context of target search kinetics [START_REF] Bénichou | Optimal reaction time for surface-mediated diffusion[END_REF][START_REF] Calandre | Accelerating search kinetics by following boundaries[END_REF][START_REF] Rupprecht | Exact mean exit time for surface-mediated diffusion[END_REF][START_REF] Monserud | Interfacial molecular searching using forager dynamics[END_REF] 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 [START_REF] Von Hippel | Facilitated target location in biological systems[END_REF][START_REF] Coppey | Kinetics of target site localization of a protein on dna: a stochastic approach[END_REF][START_REF] Berg | Diffusion-driven mechanisms of protein translocation on nucleic acids. 1. models and theory[END_REF]. It is less clear if similar optimization effects appear for transport properties. In Ref. [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF], 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 [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF] 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) [START_REF] Levesque | Taylor dispersion with adsorption and desorption[END_REF][START_REF] Berezhkovskii | Aris-taylor dispersion with drift and diffusion of particles on the tube wall[END_REF] or in fast exchange limit for non-planar geometries [START_REF] Brenner | Macrotransport processes[END_REF][START_REF] Quintard | dispersion, and interfacial transport of contaminants: Homogeneous porous media[END_REF][START_REF] Edwards | Dispersion and reaction in two-dimensional model porous media[END_REF], 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 [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF]. 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 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 u b (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 ). In our approach, we assume that the tracer particle is point-like. This assumption can be made without loss of generality since the problem of a finite size tracer particle can be studied by modifying the obstacle geometry. In this way, the particle size influences the effective diffusivity. We will also neglect inertial effects, which are irrelevant for sufficiently small tracers in viscous fluids. We can also consider a local drift field u s (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 [START_REF] Grebenkov | Imperfect diffusion-controlled reactions[END_REF], 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 [START_REF] Alexandre | Generalized taylor dispersion for translationally invariant microfluidic systems[END_REF].

The probability densities p b (r, t) and p s (r, t) to be at position r at time t in the bulk and surface obey

∂p b ∂t = ∇ • [D b ∇p b -u b (r)p b ], (1)
∂p s ∂t = ∇ s • [D s ∇ s p s -u s (r)p s ] -k d p s + k a p b , (2)
where ∇ s is the surface nabla operator. The boundary conditions, determined from probability conservation, are given by

n • [D b ∇p b -u b (r)p b ] = 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 u b (r), u s (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:

v e = x(t) -x(0) t , D e = lim t→∞ [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 Supplementary Information (SI) [54])

v e = V dr (∇x) • j b + S dS (∇ s x) • j s , (5)
j b = (u b -D b ∇)P st b , j s = (u s -D s ∇ s )P st s , (6)
where the integrals are evaluated over the volume V and the surface S of the walls in one period. The functions P st b (r) and P st s (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:

D e (F = 0) = k B T × d dF v e (F ) 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:

u b (r) = βD b F (∇x), u s (r) = βD s F (∇ s x), (8)
with β = 1/k B T . At zero force, we denote the stationary probability densities by P st,0 b , P st,0 s , both are uniform and are given by:

P st,0 b (r) = 1 V + Sδ , P st,0 s (r) = δP st,0 b , (9)
where δ = k a /k d is an adsorption length (that can be much larger than the range of interactions with the surface). The length δ quantifies how "sticky" the surface is, with surfaces becoming non-reflecting in our problem when δ is comparable to V /S, in which case the fraction of time spent on the surface become significant. 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:

P st b P st,0 b + βF f b (r), P st s P st,0 s + βF f s (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

D e =D b V dr (∇x) • ∇ P st,0 b x -f b + D s S dS (∇ s x) • ∇ s P st,0 s x -f s . (11)
The equations for f b , f s are obtained by inserting the ansatz [START_REF] Taylor | Dispersion of soluble matter in solvent flowing slowly through a tube[END_REF] into the transport equations (1, 2, 3) and expanding to linear order in F , which yields

∇ 2 f b = 0, D s ∇ 2 s [f s -P st,0 s x] = k d f s -k a f b , (12)
D b n • ∇f b = k d f s -k a f b + D b P st,0 b (∇x) • 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 [START_REF] Brenner | Macrotransport processes[END_REF], homogenization methods or Kubo formulas [START_REF] Guérin | Kubo formulas for dispersion in heterogeneous periodic nonequilibrium systems[END_REF][START_REF] Guérin | Force-induced dispersion in heterogeneous media[END_REF]. 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 → ∞ keeping k a /k d constant)
reproduce the results given in [START_REF] Brenner | Macrotransport processes[END_REF][START_REF] Quintard | dispersion, and interfacial transport of contaminants: Homogeneous porous media[END_REF][START_REF] Edwards | Dispersion and reaction in two-dimensional model porous media[END_REF]. Note that surface curvature effects arise in Eq. ( 12) via the term ∇ 2 s 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.

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 ϕ their volume fraction. An exact analytical expression for the dispersion for any ϕ is difficult to obtain even without the surface interactions. However, in the limit of small R/L (equivalently, the small ϕ limit), we can look for solutions in terms of matched asymptotic expansions, with an inner solution for f b (r) that varies at the scale |r| ∼ R , and an outer solution varying at the scale L. The auxiliary fields in the limit R → 0 are written as [START_REF] Burada | Diffusion in confined geometries[END_REF]. Squares: data obtained by numerical integration of Eqs. ( 12), [START_REF] Carusela | Computational model for membrane transporters. potential implications for cancer[END_REF]. Triangles: results of stochastic simulations of the stochastic process r(t) associated to the Fokker-Planck equations [START_REF] Marbach | Transport and dispersion across wiggling nanopores[END_REF][START_REF] Aminian | How boundaries shape chemical delivery in microfluidics[END_REF]. Circles: Monte Carlo simulation data taken from [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF], where the tracer particles are submitted to an exponentially decaying potential. The parameter δ was adjusted so that the stationary probability densities (9) in our description are the same as in the case of this potential (see SI for details).

f b (r) =      n≥0 R n+1 f (n) b (r/R, ω) (r L, inner) n≥0 R d+n F (n) b (r, ω) (r R, outer) , (14)
where r is the distance to the center of the nearest obstacle and ω 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 r L (where they are both valid). To order ϕ 4 (see SI), we find

D e = D b 1 + (dξ -1)ϕ 1 -αϕ 1 + αϕ/(d -1) + o(ϕ 4 ), α = D b + D s ξ(γ -1)(d -1) D b + D s ξ[γ(d -1) + 1] , γ = D b k a R , ξ = δ R . (15)
The above formula is exact up to the order ϕ 4 in the small ϕ limit, when ξ and γ are kept constant. It is valid for a large range of values of ϕ, as shown in Fig. 4 where it is compared to numerical solution of the full transport equations [START_REF] Bhattacharjee | Bacterial hopping and trapping in porous media[END_REF][START_REF] Carusela | Computational model for membrane transporters. potential implications for cancer[END_REF]. The theory also agrees with the Monte Carlo simulations of Ref. [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF]. This expression ( 15) is a generalization of similar Maxwell type formulas obtained for reflecting obstacles, see e.g. Refs. [START_REF] Mangeat | Effective diffusivity of Brownian particles in a two dimensional square lattice of hard disks[END_REF][START_REF] Kalnin | Calculations of the effective diffusion coefficient for inhomogeneous media[END_REF][START_REF] Maxwell | Electricity and Magnetism[END_REF][START_REF] Lebenhaft | Diffusion-controlled processes among partially absorbing stationary sinks[END_REF]. We can define a critical surface diffusivity D * s as the value of D s above which a weak attraction increases D e (so [∂D e /∂δ] δ=0 = 0 at D s = D * s ). Using Eq. ( 15), we find

D * s = D b 1 - 1 -ϕ d . (16)
Interestingly, D * s depends only on ϕ 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. [START_REF] Burada | Diffusion in confined geometries[END_REF] was found in Ref. [START_REF] Edwards | Diffusion and convective dispersion through arrays of spheres with surface adsorption, diffusion, and unequal solute partitioning[END_REF] in the fast exchange limit (γ = 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 SI). Finally, the increase in D e is rather marginal in this geometry, as can be seen in Fig. 4. 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 [START_REF] Malgaretti | Entropic transport in confined media: a challenge for computational studies in biological and soft-matter systems[END_REF][START_REF] Reguera | Entropic transport: Kinetics, scaling, and control mechanisms[END_REF]. We consider the limit of slowly varying undulations (L → ∞), where the problem can be solved using a perturbation expansion:

f b = n≥0 1 L n f (n) b x L , r ⊥ , f s = n≥0 1 L n f (n) s x L , (17)
where r ⊥ 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 SI)

D e = L→∞ 1 h d-2 (h + (d -1)δ) h 2-d D b h+(d-1)Dsδ , (18)
where

• = 1 L L 0
• 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 [START_REF] Levesque | Taylor dispersion with adsorption and desorption[END_REF]. For reflecting walls (δ = 0), we recover the well known result obtained with the Fick-Jacobs approximation, [START_REF] Jacobs | Diffusion processes[END_REF][START_REF] Zwanzig | Diffusion past an entropy barrier[END_REF]. Once again, defining D * s as the value of D s for which [∂D e /∂δ] δ=0 = 0, we find

D e (δ = 0) = D b [ h d-1 h -(d-1) ] -1
D * s = D b h d-2 h 1-d h d-1 h -d . (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. [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF].

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 δ = 0 (pure reflection) and δ → ∞ (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 (∞), then the inequality D s > D * s always holds (see SI). Hence, for any channel geometry, there is a regime for which the effective diffusivity is larger than the effective diffusivity for completely reflecting (δ = 0) or completely sticky (δ → ∞) 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. 12) and ( 13) for kaa/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 a, with a the minimal channel height), for which we obtain (see SI)

D e = L a V c V + Sδ 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 = π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. [START_REF] Martens | Hydrodynamically enforced entropic trapping of brownian particles[END_REF]. 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 [START_REF] Bénichou | Optimal reaction time for surface-mediated diffusion[END_REF][START_REF] Calandre | Accelerating search kinetics by following boundaries[END_REF][START_REF] Coppey | Kinetics of target site localization of a protein on dna: a stochastic approach[END_REF], as exemplified by the search of a gene sequence on DNA [START_REF] Von Hippel | Facilitated target location in biological systems[END_REF][START_REF] Coppey | Kinetics of target site localization of a protein on dna: a stochastic approach[END_REF][START_REF] Berg | Diffusion-driven mechanisms of protein translocation on nucleic acids. 1. models and theory[END_REF]. 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. [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF] 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. [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF] 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.

Computer time for this study was provided by the computing facilities MCIA (Mesocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour.

SUPPLEMENTAL MATERIAL

Here, we provide details of calculations that support the results of the main text:

• a proof of the general expression of the average velocity v e (Section A),

• some reminders of tensor analysis for diffusion on surfaces (Section B),

• a description of stochastic simulations of the process r(t) (Section C)

• a derivation of the expression of D e in periodic arrays of spherical obstacles, including a comparison with the results of the literature (Section D),

• a derivation of the dispersion for slowly varying and fast varying channel widths (Sections E and F).

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

x(t) -x 0 = Ω dr(x -x 0 )p b (r, t|r 0 , i 0 ) + Σ dS(r)(x -x 0 )p s (r, t|r 0 , i 0 ) , (A1)
where Ω and Σ are the volume and surface of the whole medium. Taking the temporal derivative and using Eqs. (1,2) yields

∂ t x(t) = Ω dr (x -x 0 )∇ • [(D b ∇ -u b )p b (r, t|r 0 , i 0 )] + Σ dS(r)(x -x 0 ) {∇ s • [(D s ∇ s -u s )p s (x, t|r 0 , i 0 )] -k d p s (r, t|r 0 , i 0 ) + k a p b (r, t|r 0 , i 0 )} . (A2)
After performing integrations by parts and using the boundary condition Eq. ( 3), we obtain

∂ t x(t) = Ω dx(∇x) • (u b -D b ∇) p b (r, t|r 0 , i 0 ) + Σ dS(x)(∇ s x) • (u s -D s ∇ s ) p s (r, t|r 0 , i 0 ) . (A3)
We now introduce P s (r, t) and P b (r, t), which are the probability densities at t to be respectively on the surface or in the bulk, at position r modulo the period, i.e.:

P b (r, t) = L p b (r + L, t) , P s (r, t) = L p s (r + L, t) , (A4)
where the vectors L are all the vectors such that the environments at r and r + L are identical (for example, in the channel geometry L = nLe x with n integer). In Eq. (A3), ∇x, ∇ s x, u b and u s are periodic in r. Therefore, splitting the integral into multiple integrals (each of them covering one period) and using the above relation yields

∂ t x(t) = V dr (∇x) • (u b -D b ∇) P b (r, t|r 0 , i 0 ) + S dS(r) (∇ s x) • (u s -D s ∇ s ) P s (r, t|r 0 , i 0 ) , (A5)
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 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 st i (r) = V dr 0 P i (r, t|r 0 , b)P st b (r 0 ) + S dS(r 0 )P i (r, t|r 0 , s)P st s (r 0 ) ∀t , (i ∈ {b, s}) (A6)
so that averaging Eq. (A5) over initial conditions and integrating over t leads to

v e = ∂ t x(t) = x(t) -x(0) t = V dr (∇x) • (u b -D b ∇) P st b (r) + S dS(r) (∇ s x) • (u s -D s ∇ s ) P st s (r) , (A7)
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. (A5).

Appendix B: Tensor analysis for surface diffusion

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

g ij = e s,i • e s,j , e s,i = ∂r s ∂y i , (B1)
where the suffix s indicates that the position remains on the surface. We use the notation |g| = det(g ij ) and g ij for the inverse of the metric tensor (i.e. g ij g jk = δ 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 φ, which are well known and are given by

dS = dy 1 dy 2 ...dy d-1 |g| , ∇ s φ = g ij ∂φ ∂y j e s,i , ∇ 2 s φ = 1 |g| ∂ ∂y i |g| g ij ∂φ ∂y j . ( B2
)
Appendix C: Algorithm for stochastic simulations

We have used the following algorithm to simulate the stochastic process r(t) associated with 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 = √ 2D b dW i (with dW i = 0 and dW i dW j = dtδ ij ) during the time step dt in each direction i = {x, y, z}.

• At the end of a step, if the tracer particle ends up inside the obstacle, it switches to the attached state with probability P adsorb = k a πdt/D b and it is reflected with probability 1 -P adsorb , see Ref. [START_REF] Singer | Partially reflected diffusion[END_REF].

• An attached particle switches to the detached state with probability 1 -e -k d dt .

• If the particle is attached to a spherical 3D obstacle, the stochastic equation (with Ito convention) for the angular coordinates (θ, φ) is

dθ t = D s dt R 2 tan θ t + √ 2D s R dW θ , dφ t = 2D s R 2 sin 2 θ t dW φ , dW 2 φ = dW 2 θ = dt. (C1)
This equation can be read off from the Fokker-Planck equation for p s (r, t), which in terms of spherical angles (θ, φ) reads (in absence of detachment event) :

∂ t p s (θ, φ, t) = D s R 2 sin θ ∂ θ [sin θ∂ θ p s ] + D s R 2 sin 2 θ ∂ 2 φ p s . (C2)
Here, p s dS = p s R 2 sin θdθdφ is the probability to find the particle with the angular coordinates within the interval [θ, θ + dθ] × [φ, φ + dφ]. We now consider q(θ, φ) the probability distribution function of (θ, φ) so that q(θ, φ)dθdφ is the probability to observe the particle in [θ, θ + dθ] × [φ, φ + dφ]. We thus have q = p s R 2 sin θ and q obeys

∂q ∂t = ∂ 2 ∂φ 2 D s R 2 sin 2 θ q + ∂ ∂θ ∂ ∂θ D s R 2 q - D s R 2 tan θ q . (C3)
It is clear that the stochastic differential equation that corresponds to this Fokker-Planck equation is Eq. (C1).

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

2D obstacles

Here, we compute D e for a square array of 2D spheres of radius R. The volume fraction is given by ϕ = πR 2 /L 2 . We chose the units of length so that L = 1, without loss of generality. We use polar coordinates (r, θ), with r the distance to the center of the disk and θ the angle with the x-axis. The surface is parametrized by the coordinate y 1 = θ. Using the formulas of Section B, we obtain ∇ 2 s x = -cos θ/R. Equations [START_REF] Bhattacharjee | Bacterial hopping and trapping in porous media[END_REF][START_REF] Carusela | Computational model for membrane transporters. potential implications for cancer[END_REF] written in polar coordinates are: . After an integration by parts in Eq. ( 11), we find that the effective diffusion constant is

               (bulk) 1 r ∂ ∂r r ∂f b ∂r + 1 r 2 ∂ 2 f b ∂θ 2 = 0 , (r > R) (surface) D s R 2 ∂ 2 f s ∂θ 2 = k d f s (θ) -k a f b (r = R, θ) -D s P st,0 s cos θ R , (r = R) (boundary condition) -D b ∂f b ∂r (r = R, θ) = k d f s -k a f b -D b P st,0 b cos θ . (r = R) (
D e = D b (1 -ϕ)P st,0 b + D s ϕξP st,0 b + D b R 2π 0 dθ cos(θ)f b (r = R, θ) -D s 2π 0 dθ cos(θ)f s (θ) . ( D2
)
We now introduce the dimensionless parameters ξ = δ/R and γ = D b /(k a R). We wish to compute D e in the limit R → 0 (or, equivalently, the limit of vanishing volume fraction ϕ → 0), when ξ and γ are kept constant. We solve this problem by using strongly localized perturbation analysis [START_REF] Ward | Strong localized perturbations of eigenvalue problems[END_REF], considering that the solution f b varies at two different length scales r ∼ R and r ∼ 1 when R → 0 (we recall that the period is L = 1). We introduce an inner solution (valid near the obstacles), depending on r = r/R and θ, which has an expansion in the limit R → 0 at fixed r = r/R:

f b (r) = P st,0 b R f 0 b (r, θ) + R 2 f 1 b (r, θ) + R 4 f 2 b (r, θ) + ... , f s (r) = P st,0 b R 2 f 0 s (θ) + R 2 f 1 s (θ) + R 4 f 2 s (θ) + ... . (D3)
Note that we have introduced the factor P st,0 b since it is obvious in Eq. (D1) that the solutions are proportional to this factor. We introduce an outer solution, far from the obstacle, which admits the following expansion when R → 0 at fixed r:

f b (r) = P st,0 b R 2 F 0 b (r, θ) + R 2 F 1 b (r, θ) + R 4 F 2 b (r, θ) + ... . (D4)
It follows that all f i b (r, θ) and F i b (r, θ) are harmonic functions for all orders i, and that all F i b (r, θ) must satisfy periodic conditions on the sides of the unit square. The equations for f i s and for the value of f i b near the obstacles are found by inserting Eq. (D3) into Eq. (D1) and identifying terms order by order, leading to:

D s γξ ∂ 2 f i s ∂θ 2 = D b f i s -ξf i b -δ i,0 ξ 2 D s γ cos θ, -ξγ ∂f i b ∂ r (r = 1, θ) = f i s -ξf i b -δ i,0 ξγ cos θ, (D5)
where δ i,0 is the Kronecker delta symbol. Additional conditions to determine f i b and F i b come from the requirement that the inner and the outer solution coincide in their regime of common validity R r 1, so that the expansion of Eq. (D3) for r → ∞ coincides with the expansion of (D4) for r → 0.

At leading order, f 0 b is the harmonic function that does not diverge for large r and that satisfies Eq. (D5) for i = 0:

f 0 b (r, θ) = A 0 r cos(θ) must match F 0 b , f 0 s (θ) = f 0 s cos(θ), A 0 = -α, f 0 s = ξ D s ξ(γ + 1) -D b D s ξ(γ + 1) + D b , α = D b + D s ξ(γ -1) D b + D s ξ(γ + 1) . (D6)
We see that f 0 b ≈ A 0 cos θ/r for large r, so that f b A 0 cos θR 2 /r in the region R r 1. This implies that the leading order of the outer solution must be of order R 2 , as assumed in Eq. (D4), and that F 0 b (r → 0, θ) = A 0 cos θ/r, as indicated in Eq. (D6). Since F 0 b is also harmonic with periodic boundary conditions, it can be expressed as

F 0 b (r, θ) = -2πA 0 e x • ∇G , ( D7
)
where G is the pseudo Green's function for the unit-square with periodic conditions, which satisfies the equation:

-∇ 2 G = δ(x)δ(y) -1 , G(x ± 1, y) = G(x, y ± 1) = G(x, y) . ( D8
)
The Green's function G is given in closed form in Ref. [START_REF] Mamode | Fundamental solution of the laplacian on flat tori and boundary value problems for the planar poisson equation in rectangles[END_REF], its behaviour near the origin reads

G(r, θ) = - log r 2π + r 2 4 + C - C 0 8π r 4 cos(4θ) + O(r 8 ) , (D9)
with C 0 a constant that we will not need to detail here. Eq. (D7) is justified by noting that F 0 b near the origin is

F 0 b (r, θ) = A 0 r cos(θ) matching f 0 b -πA 0 cos(θ)r must match f 1 b + C 0 cos(3θ)r 3 must match f 2 b + O(r 7 ) , (D10)
where we see that the 1/r term correctly matches the corresponding 1/r term in Eq. (D6). We have also indicated that the terms of order r and r 3 will have to be matched by corresponding terms r and r3 in the expansions of f 1 b and f 2 b , respectively. In particular, f 1 b (r, θ) -πA 0 r cos θ for large r, since f 1 b is harmonic, it is of the form

f 1 b (r, θ) = A 1 r cos(θ) must match F 1 b + B 1 r cos(θ) matching F 0 b , B 1 = -πA 0 . (D11)
Solving Eq. (D5) for i = 1, we obtain f 1 s and the constant A 1 :

A 1 = αB 1 , f 1 s (θ) = f 10 s cos(θ) , f 10 s = 2ξD b B 1 D s ξ(γ + 1) + D b . (D12)
The outer solution at next-to-leading order F 1 b and its expansion near the origin then reads:

F 1 b (r, θ) = -2πA 1 e x • ∇G = A 1 r cos(θ) matching f 1 b -πA 1 cos(θ)r must match f 2 b + C 1 cos(3θ)r 3 must match f 3 b +O(r 7 ) . ( D13
)
Second order: We continue the iteration and we look for f 2 b in the form

f 2 b (r, θ) = A 2 r cos(θ) must match F 2 b + B 2 r cos(θ) matching F 1 b + K 2 r3 cos(3θ) must match F 3 b + C 0 r3 cos(3θ) matching F 0 b , B 2 = -πA 1 , (D14)
so that it matches corresponding terms for F 1 b and F 0 b in Eqs. (D10) and (D13), respectively. The expression for A 2 and f 2 s can be found by solving Eq. (D5) with i = 2:

A 2 = αB 2 , f 2 s (θ) = f 20 s cos(θ) + f 21 s cos(3θ) , f 20 s = 2ξD b B 2 D s ξ(γ + 1) + D b . (D15)
The values of K 2 and f 21 s are not written because they do not contribute to D e . For the outer solution, we obtain

F 2 b (r, θ) = -2πA 2 e x • ∇G = A 2 r cos(θ) matching f 2 b -πA 2 cos(θ)r must match f 3 b + C 2 cos(3θ)r 3 must match f 4 b + O(r 7 ) . (D16)
Third order: At this order, the expression of the inner solution reads

f 3 b (r, θ) = A 3 r cos(θ) must match F 3 b + B 3 r cos(θ) matching F 2 b + K 3 r3 cos(3θ) must match F 4 b + C 1 r3 cos(3θ) matching F 1 b , B 3 = -πA 2 . ( D17
)
This form is chosen so that the corresponding terms in F 2 b and F 1 b are matched. Identically to Eq. (D15), we get

A 3 = αB 3 , f 3 s (θ) = f 30 s cos(θ) + f 31 s cos(3θ) , f 30 s = 2ξD b B 3 D s ξ(γ + 1) + D b . (D18)
At this order of perturbation, we find that the effective diffusivity as given in Eq. (D2) is given by:

D e = P st,0 b D b (1 -ϕ) + D s ξϕ + D b ϕ A 0 + 3 i=1 R 2i (A i + B i ) -D s ϕ 3 i=0 R 2i f i s + o(R 8 ) . (D19)
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 st,0 b 1 -2αϕ + 2α 2 ϕ 2 -2α 3 ϕ 3 + 2α 4 ϕ 4 + o(ϕ 4 ) = D b 1 + (2ξ -1)ϕ 1 -αϕ 1 + αϕ + o(ϕ 4 ) . (D20)

3D spherical obstacles

We now consider a 3-dimensional cubic array of spheres of radius R. Here, we use spherical coordinates (r, θ, φ), 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 (θ = π/2, φ = nπ/2), and θ = 0 and θ = π. The long-time diffusion constant is given by

D e = P st,0 b [D b (1 -ϕ) + 2ϕξD s ] + π 0 dθ 2π 0 dφ sin(θ) cos(θ) D b R 2 f b -2D s Rf s , (D21)
where ϕ = 4π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:

f b (r) = P st,0 b R f 0 b (r, θ, φ) + Rf 1 b (r, θ, φ) + R 2 f 2 b (r, θ, φ) + ... , f s (r) = P st,0 b R 2 f 0 s (θ, φ) + Rf 1 s (θ, φ) + R 2 f 2 s (θ, φ) + ... , (R → 0, fixed r = r/L, θ, φ) (D22)
where we have factorized by

P st,0 b = [1 + (3ξ -1)ϕ] -1 . The outer solution is f b (r) = P st,0 b R 3 F 0 b (r, θ, φ) + RF 1 b (r, θ, φ) + R 2 F 2 b (r, θ, φ) + ... . (R → 0, fixed r, θ, φ) (D23)
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 i b and F i b are harmonic functions, (ii) F i b satisfy periodic boundary conditions, (iii) the behavior of the inner and the outer solution is the same at all orders in the region R r 1, and (iv) the inner functions satisfy

       (surface) γD s ξ 1 sin(θ) ∂ ∂θ sin(θ) ∂f i s ∂θ + 1 sin(θ) 2 ∂ 2 f i s ∂φ 2 = D b f i s -ξf i b -2δ i,0 ξ 2 D s γ cos(θ) , (boundary condition) -ξγ ∂f i b ∂ r (r = 1, θ) = f i s -ξf i b -δ i,0 ξ cos(θ) . ( D24
)
At leading order, the inner solution is determined from Eq. (D24) and reads

f 0 b (r, θ) = A 0 r2 cos(θ) must match F 0 b , f 0 s (θ) = f 0 s cos(θ), A 0 = - α 2 , α = 2D s ξ(γ -1) + D b D s ξ(2γ + 1) + D b , f 0 s = ξ[2D s ξ(2γ + 1) -D b ] 2[D s ξ(2γ + 1) + D b ] . (D25)
It is useful to introduce the pseudo-Green's function G with periodic boundary conditions, which satisfies

-∇ 2 G = δ(r) -1 , G(x ± 1, y, z) = G(x, y ± 1, z) = G(x, y, z ± 1) = G(x, y, z). (D26)
In fact, we will not need the full expression of G but only its behavior near the origin, namely

G(r, θ, φ) = r→0 1 4πr + r 2 6 + r 4 A 4 (θ, φ) + r 8 A 8 (θ, φ) + ... ( D27
)
with A i (θ, φ) functions that are built from Legendre polynomials and respect the symmetry of the problem. For example, A 4 is built as a combination aP 0 4 (cos θ) + bP 4 4 (cos θ) cos(4φ) (with P m n the Legendre polynomials) with the requirement that, for φ = nπ/2 (with n an integer), A 4 does not depend on cos(2θ). The only possible form for A 4 is

A 4 (θ, φ) = C 0 {20 cos(2θ) + 35 cos(4θ) + 9 -5 cos(4φ)[4 cos(2θ) -cos(4θ) -3]} (D28)
and C 0 a constant that we will not need to determine here. We claim that F 0 b (r, θ) = -4πA 0 e z • ∇G. Indeed, this is a harmonic function with periodic boundary conditions, and its expansion near the origin reads:

F 0 b (r, θ) = -4πA 0 e z • ∇G = A 0 r 2 cos(θ) matching f 0 b - 4π 3 A 0 cos(θ)r must match f 3 b + D 0 [3 cos(θ) + 5 cos(3θ)]r 3 must match f 5 b + E 0 (θ, φ)r 7 must match f 9 b + o(r 7 ) , (D29)
with D 0 = -128πC 0 and E 0 (θ, φ) a combination of Legendre polynomials, for which we will only need the property:

π 0 dθ 2π 0 dφ sin(θ) cos(θ)E 0 (θ, φ) = 0 . (D30)
Next, at first and second order, we obtain the trivial solutions

f 1 b = f 1 s = F 1 b = 0 and f 2 b = f 2 s = F 2 b = 0.
The inner solution at third order is obtained by solving Eq. (D24) with i = 3 and requiring that the large distance behavior matches with the corresponding term in the expansion of F 0 b in Eq. (D29), we obtain

f 3 b (r, θ) = A 3 r2 cos θ must match F 3 b + B 3 r cos θ matching F 0 b , f 3 s (θ) = cos θ f 3 s , B 3 = - 4π 3 A 0 , A 3 = B 3 2 α, f 3 s = 3ξD b B 3 2 (D s ξ(2γ + 1) + D b ) . (D31)
The outer solution at this order and its expansion near the origin is 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 5 b = 0. We continue the iteration and we find at order 6:

F 3 b (r, θ, φ) = -4πA 3 e z • ∇G = A 3 r 2 cos(θ) matching f 3 b - 4π 3 A 3 cos(θ)r must match f 6 b + D 1 [3 cos(θ) + 5 cos(3θ)]r 3 must match f 8 b + O(r 7 ) . (
f 6 b = A 6 r2 cos(θ) must match F 6 b + B 6 r cos(θ) matching F 3 b , f 6 s (θ) = cos(θ) f 6 s , B 6 = - 4π 3 A 3 , A 6 = B 6 2 α , f 6 s = 3ξD b B 6 2 (D s ξ(2γ + 1) + D b )
, (D35)

F 6 b (r, θ, φ) = -4πA 6 e z • ∇G = A 6 r 2 cos(θ) matching f 6 b - 4π 3 A 6 cos(θ)r must match f 9 b + D 2 [3 cos(θ) + 5 cos(3θ)]r 3 must match f 11 b + O(r 7 ) . (D36)
For the inner solution at 7 th order, there is no term to be matched, so f 7 b = f 7 s = 0. However, for the outer solution F 7 b , we have to match a term from f 5 b . We can build F 7 b from (e z • ∇) 2 F 0 b so that its expansion near the origin reads

F 7 b (r, θ) = 4A 5 3A 0 (e z • ∇) 2 F 0 b = A 5 r 4 [3 cos(θ) + 5 cos(3θ)] matching f 5 b + 64A 5 D 0 A 0 r cos(θ) must match f 10 b + O(r 3 ) . (D37)
8 th order: Here, the form of the inner solution is

f 8 b (r, θ) = [3 cos(θ) + 5 cos(3θ)] A 8 r4 must match F 10 b + [3 cos(θ) + 5 cos(3θ)] D 1 r3 matching F 3 b , f 8 s (θ) = f 8 s [3 cos(θ) + 5 cos(3θ)] . (D38)
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 8 b = 0. Finally, the inner solution at 9 th order reads

f 9 b (r, θ) = A 9 r2 cos(θ) must match F 9 b + B 9 r cos(θ) matching F 6 b + E 0 (θ, φ)r 7 matching F 0 b + C 9 E 0 (θ, φ) r8 must match F 15 b , (D39)
f 9 s (θ) = f 9 s cos(θ) , A 9 = B 9 2 α , B 9 = - 4π 3 A 6 , f 9 s = 3ξD b B 9 2 (D s ξ(2γ + 1) + D b ) . ( D40
)
Since E 0 satisfies Eq. (D30), the two last terms in f 9 b do not contribute in the expression of D e . We stop here the iteration procedure and we consider the long-time dispersion Eq. (D21) which reads:

D e = P st,0 b {D b (1 -ϕ) + 2D s ξϕ + D b ϕ A 0 + 3 n=1 R 3n (A 3n + B 3n ) -2D s ϕ 3 n=0 R 3n f 3n s + o(R 12 ) . (D41)
Using the above results for A n , B n , f n s , this expression can be simplified, leading to a re-summation and to the result

D e = D b P st,0 b 1 - 3ϕα 2 1 - αϕ 2 + αϕ 2 2 - αϕ 2 3 + o(ϕ 4 ) = D b 1 + (3ξ -1)ϕ 1 -αϕ 1 + αϕ 2 + o(ϕ 4 ) .
(D42)

Comparison with existing results

a. Monte Carlo data of Ref. [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF] Here, we compare the prediction of Eq. ( 15) with Monte Carlo simulations of Ref. [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF], as they appear in Figure 9 of their Supplementary Information, 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. [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF], the tracer particle is subject to an exponentially decaying potential U (D), with D the distance to the crowder's wall. The link between the parameter δ in our description and the potential U (D) can be made by considering the simple one-dimensional situation in which the tracer particle diffuses between x = 0 and x = H and is subject to the potential U (x), assumed to vanish for distances larger than a length λ, with λ H. In this case, the probability p Σ to find the particle at a distance from the surface less than ε (with λ H) reads

p Σ = 0 e -βU (x) dx H 0 e -βU (x) dx = 0 [e -βU (x) -1]dx + ε H 0 [e -βU (x) -1]dx + H ∞ 0 e -βU (x) -1 dx ∞ 0 e -βU (x) -1 dx + H . (D43)
In our formalism, we readily find p Σ = δ/(δ + H), so that the expression of the length δ in terms of U is

δ = ∞ 0 e -βU (D) -1 dD . (D44)
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 → ∞, leading to γ = 0 in Eq. ( 15). Finally, Ref. [START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF] considers a tracer particle of finite radius r d in presence of spherical crowders of radius r c and volume fraction φ, 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 ϕ = φ(1 + r d /r c ) 3 . We obtain ϕ = 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.

b. Comparison with the result of macrotransport theory

In Ref. [START_REF] Edwards | Diffusion and convective dispersion through arrays of spheres with surface adsorption, diffusion, and unequal solute partitioning[END_REF], 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 (r s , t) = k ed p b (r s , t) and p inside (r s , t) = K ed p b (r s , t) where k ed and K ed are Henry's laws constants, and p inside (r, t) is the probability density of tracer position inside the spheres. In principle, if one takes K ed = 0 to prevent the tracers from entering into the obstacles, the physical situation of Ref. [START_REF] Edwards | Diffusion and convective dispersion through arrays of spheres with surface adsorption, diffusion, and unequal solute partitioning[END_REF] is identical to our formalism in the fast exchange limit (γ = 0). The result of Ref. [START_REF] Edwards | Diffusion and convective dispersion through arrays of spheres with surface adsorption, diffusion, and unequal solute partitioning[END_REF] for a periodic array of spheres in 3D is given by their equation (4.23) and reads with our notations:

D e = D b 1 + (3ξ -1)ϕ   1 - 3ϕ 2 1 + Ds D b ξ + ϕ 1 -2 Ds D b ξ   + o(ϕ 4 ) = D b 1 + (3ξ -1)ϕ 1 -ϕ 1 + αϕ 2 + o(ϕ 4 ) , (D45)
which is identified by writing K ed = γ ed = 0, k ed α ed = 3ξϕ, Γ ed = D s ξ/D b (the suffix "ed" refers to notations of Ref. [START_REF] Edwards | Diffusion and convective dispersion through arrays of spheres with surface adsorption, diffusion, and unequal solute partitioning[END_REF]). However, the above equation (D45) disagrees with our Eq. ( 15), since the numerator contains a term 1 -ϕ rather than 1 -αϕ.

In Fig. 4, we show our prediction [Eq. [START_REF] Burada | Diffusion in confined geometries[END_REF]] along with the result (D45), compared with the numerical solution of our formalism. Our analytical formula agrees with the numerical results, while the formula D45 does not and in particular the non-monotonicity of D e with δ is absent. An additional check of the validity of our formalism is provided by comparison with direct stochastic numerical simulations of the process r(t) corresponding to the Fokker-Planck equations ( 1) and (2) (described in Section 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. Furthermore, the above results for the diffusivity Eqs. (D20) and (D42) are in surprisingly good agreement with simulations even for higher values of ϕ, as shown in Fig. 5 in the case of spherical obstacles.

13) read

                     ∂ 2 f b ∂x 2 + ∂ 2 f b ∂y 2 = 0 , [0 < y < h(x)] D s 1 + h (x) 2 ∂ ∂x 1 1 + h (x) 2 ∂f s ∂x -k d f s + k a f b = -D s P st,0 s h (x)h (x) [1 + h (x) 2 ] 2 , [y = h(x)] D b 1 + h (x) 2 ∂f b ∂y -h (x) ∂f b ∂x = k d f s -k a f b - D b P st,0 b h (x) 1 + h (x) 2 . [y = h(x)] (E1)
Furthermore, f b and f s are periodic with period L, and satisfy ∂ y f b | y=0 = 0 (from the symmetry y → -y). The stationary probability densities and the effective diffusivity are given by

P st,0 b = 1 V + Sδ , P st,0 s = δP st,0 b , V = 2 L 0 dx h(x) , S = 2 L 0 dx 1 + h (x) 2 , ( E2
)

D e = P st,0 b V D b + 2D s P st,0 s L 0 dx 1 + h (x) 2 -D b L 0 dx h(x) -h(x) dy ∂f b ∂x -2D s L 0 dx 1 + h (x) 2 ∂f s ∂x . ( E3
)
We now analyze the above equations in the limit L → ∞, when all other parameters are held constant. Here, we write h(x) = ζ(x/L), where the function ζ(x) does not depend on L, and we use the notation x = x/L. The expansion of f b takes the form:

f b (r) = f 0 b (x, y) + L -1 f 1 b (x, y) + L -2 f 2 b (x, y) + ... , f s (x) = f 0 s (x) + L -1 f 1 s (x) + L -2 f 2 s (x, y) + ... ( E4
)
Note that this expansion in powers of L is equivalent to an expansion in terms of the small parameter ε = 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 x = x/L:

                     1 L 2 ∂ 2 f b ∂ x2 + ∂ 2 f b ∂y 2 = 0 , [0 < y < ζ(x)] D s 1 + ζ (x) 2 /L 2 1 L 2 ∂ 2 fs ∂ x2 - 1 L 4 ζ (x)ζ (x) 1 + ζ (x) 2 /L 2 ∂ fs ∂ x = k d f s -k a f b - D s P st,0 s ζ (x)ζ (x) L 3 (1 + ζ (x) 2 /L 2 ) 2 , [y = ζ(x)] D b 1 + ζ (x) 2 /L 2 ∂f b ∂y - ζ (x) L 2 ∂f b ∂ x [x, h(x)] = k d f s -k a f b - D b P st,0 b ζ (x) L 1 + ζ (x) 2 /L 2 . [y = ζ(x)] (E5) P st,0 b = 1 CL 1 + O L -2 , P st,0 s = δ CL 1 + O L -2 , C = 2 ( ζ + δ) . (E6)
We insert the ansatz (E4) into the above equations to find successive orders 1/L n . At leading order, we obtain

∂ 2 f 0 b ∂y 2 = 0 , (f 0 s -δf 0 b ) y=ζ(x) = 0 , D b k d ∂f 0 b ∂y -f 0 s + δf 0 b y=ζ(x) = 0 . (E7)
Taking into account the condition ∂ y f 0 b | y=0 = 0, the above equations indicate that f 0 b and f 0 s do not depend on y, and that they are proportional at the surface:

f 0 b (x, y) = f 0 b (x) , f 0 s (x) = δf 0 b (x) . (E8)
Next, at first order in 1/L, we obtain the same equations as (E8), leading to

f 1 b (x, y) = f 1 b (x) and f 1 s (x) = δf 1 b (x). At second order, we obtain ∂ 2 f 0 b ∂ x2 + ∂ 2 f 2 b ∂y 2 = 0 , [0 < y < ζ(x)] (E9) D s ∂ 2 f 0 s ∂ x2 -k d f 2 s + k a f 2 b = 0 , [y = ζ(x)] (E10) D b ∂f 2 b ∂y - ζ (x) 2 2 ∂f 0 b ∂y -ζ (x) ∂f 0 b ∂ x = k d f 2 s -k a f 2 b - D b C ζ (x) . [y = ζ(x)] (E11)
We see that Eq. (E9) simplifies to

∂f 2 b ∂y (x, y) = -f 0 b (x)y , (E12)
where we have used the symmetry y → -y. Next, combining Eqs. (E10), (E11) and (E12) and using the result at order 0 [Eq. (E8)], we obtain

f 0 b (x) (D b ζ(x) + D s δ) + f 0 b (x)D b ζ (x) = D b ζ (x) C . ( E13
)
The solution of this first-order differential equation that satisfies the periodic boundary conditions is

f 0 b (x) = 1 C (D b ζ(x) + D s δ)   D b ζ(x) - D b ζ/ (D b ζ + D s δ) (D b ζ + D s δ) -1   . (E14)
Inserting this result into Eq. (E3), we find (after a few lines of algebra)

D e = 1 h + δ (D b h + D s δ) -1 + O(L -2 ) . ( E15
)
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 θ the angle around this axis. due to rotational invariance, f b (r, x, θ) and f s (x, θ) do not depend on θ. We parametrize the surface with x and θ. The metric is identified by calculating e s,x = ∂ x r s = e x + h (x)e r and e s,θ = h(x)e θ with e r the unit vector perpendicular to the channel axis, and e θ = e x × e r . We thus obtain g xx = 1 + h (x) 2 , g x,θ = 0 and g θθ = h(x) 2 . Using Section B, Eqs. [START_REF] Bhattacharjee | Bacterial hopping and trapping in porous media[END_REF][START_REF] Carusela | Computational model for membrane transporters. potential implications for cancer[END_REF] become

                   1 r ∂ ∂r r ∂f b ∂r + ∂ 2 f b ∂x 2 = 0 , D s 1 + h (x) 2 ∂ 2 f s ∂x 2 + h (x) h(x) - h (x)h (x) 1 + h (x) 2 ∂f s ∂x = k d f s -k a f b + D s P st,0 s 1 + h (x) 2 h (x) h(x) - h (x)h (x) 1 + h (x) 2 , D b 1 + h (x) 2 ∂f b ∂r -h (x) ∂f b ∂x [x, h(x)] = k d f s -k a f b - D b P st,0 b h (x) 1 + h (x) 2 , (E16)
and D e is given by

D e = P st,0 b V D b + 2πD s P st,0 s L 0 dxh(x) 1 + h (x) 2 -2πD b L 0 dx h(x) 0 dyy ∂f b ∂x -2πD s L 0 dxh(x) 1 + h (x) 2 ∂f s ∂x . (E17)
As in the 2D case, we set h(x) = ζ(x/L) and we use the notation x = x/L, which gives

                   1 r ∂ ∂r r ∂f b ∂r + 1 L 2 ∂ 2 f b ∂ x2 = 0 , D s L 2 ζ (x) 2 ∂ 2 f s ∂ x2 + ζ (x) ζ(x) - ζ (x)ζ (x) L 2 + ζ (x) 2 ∂f s ∂ x = k d f s -k a f b + LD s P st,0 s L 2 + ζ (x) 2 ζ (x) ζ(x) - ζ (x)ζ (x) [L 2 + ζ (x) 2 ] , D b 1 + ζ (x) 2 /L 2 ∂f b ∂r - ζ (x) L 2 ∂f b ∂ x [x,ζ(x)] = k d f s -k a f b - D b P st,0 b ζ (x) L 1 + ζ (x) 2 /L 2 , (E18)
with

P st,0 b = 1 C L 1 + O L -2 , P st,0 s = δ C L 1 + O L -2 , C = π h 2 + 2πδ h . (E19)
In the limit L → ∞ 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 0 b (x, r) = f 0 b (x) , f 0 s (x) = δf 0 b (x) , f 1 b (x, r) = f 1 b (x) , f 1 s (x) = δf 1 b (x) . (E20)
The equations at 2 nd order are

1 r ∂ ∂r r ∂f 2 b ∂r + ∂ 2 f 0 b ∂ x2 = 0 , [0 < r < ζ(x)] (E21) D s ∂ 2 f 0 s ∂ x2 + ζ (x) ζ(x) ∂f 0 s ∂ x = k d f 2 s -k a f 2 b + D s δζ (x) C ζ(x) , [r = ζ(x)] (E22) D b ∂f 2 b ∂r - ζ (x) 2 2 ∂f 0 b ∂r -ζ (x) ∂f 0 b ∂ x = k d f 2 s -k a f 2 b - D b ζ (x) C . [r = ζ(x)] (E23)
The integration of Eq. (E21) leads to ∂ r f

2 b = - rf 0 b (x) 2
which can be combined with Eqs. (E22, E23), leading to

f 0 b (x) D b ζ 2 (x) 2 + D s δζ(x) + f 0 b (x) [D b ζ(x)ζ (x) + D s δζ (x)] = 1 C [D b ζ(x)ζ (x) + D s δζ (x)] . (E24)
Solving this differential equation, we obtain:

f 0 b (x) = 1 C   1 - 1 ζ(x) (D b ζ(x) + 2D s δ) h -1 (D b h + 2D s δ) -1   . (E25)
Inserting this result into Eq. (E17) yields the D e at leading order in 1/L:

D e = 1 h (h + 2δ) h -1 (D b h + 2D s δ) -1 + O L -2 .
(E26)

Upper bounds on critical surface diffusion D * s

Here we show that

D * s ≤ D 0 s ≤ D b , (E27)
where we recall that D * s is the value of D s for which ∂D e /∂δ = 0, and D 0 s is the value of D s for which D e (δ = 0) = D e (δ = ∞). Using Eq. ( 18) it is straight forward to show that

D * s = D b h d-2 h 1-d h d-1 h -d , D 0 s = D b h d-2 h 2-d h d-1 h 1-d . (E28)
Here, we apply the Cauchy-Schwarz inequality which states that, for all integrable functions f, g:

f g 2 ≤ f 2 g 2 , ( E29
)
The choice f = h 1-d/2 and g = h -d/2 leads to the inequality h 1-d 2 ≤ h -d h 2-d , which leads to D * s ≤ D 0 s . Now, using again Eq. (E29), we find h d-2 2 ≤ h d-1 h d-3 and h 2-d 2 ≤ h 1-d h 3-d , leading to

h d-2 h 2-d h d-1 h 1-d ≤ h d-3 h 3-d h d-2 h 2-d . ( E30
)
For d = 2, the inequality 1 ≤ h h -1 directly yields D 0 s ≤ D b . For higher dimensions, we find by induction

D 0 s D b = h d-2 h 2-d h d-1 h 1-d ≤ h d-3 h 3-d h d-2 h 2-d ≤ ... ≤ 1 h h -1 ≤ 1 , (E31)
so we recover Eq. (E27).

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 where there is a uniform potential v s . In terms of the coordinates perpendicular to the channel axis, the overall potential V (r ⊥ ) 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) = C(x)
dr ⊥ e -βV (r ⊥ ) = s d-1 h(x) d-1 + a d-1 h(x) d-2 exp(-βv s ) , (E32)

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 → 0 with exp(-βv s ) = δ [see Eq. which is precisely Eq. ( 18) of our Letter.

Appendix F: Dispersion in highly undulated channels (a L) -Proof of Eq. ( 20)

FIG. 6. Schematic drawing of a channel where the half-width w is shown.

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. 6. For simplicity, we assume that the corrugation is symmetric about x = 0. This assumption can in fact 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 e s,r = e r + w (r)e x which leads to g rr = 1 + w 2 (r). For an axisymmetric 3D channel, we parametrize the surface with r and θ, the components of the metric read g rr = 1 + w 2 (r), g rθ = 0 and g θθ = r 2 . Here, we consider the large corrugation limit a → ∞, by assuming that w(r) = w(r/a), where w does not depend on a. In terms of the rescaled variable r = r/a, equations [START_REF] Bhattacharjee | Bacterial hopping and trapping in porous media[END_REF][START_REF] Carusela | Computational model for membrane transporters. potential implications for cancer[END_REF] read (for d = 2 or d = 3): where we used the notation fb = f b /P st b,0 and fs = f s /P st b,0 . At leading order in the limit a → ∞, we find with V p the volume of the peripheral region. Using P st b,0 = [V + δS] -1 , we thus recover Eq. ( 20).

                           (bulk)
∂ 2 f 0 b ∂x 2 = 0 , k d f 0 s -k a f 0 b [ w(r

FIG. 1 .

 1 FIG. 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).

FIG. 2 .

 2 FIG. 2. (color online) Long-time dispersion for a periodic array of spherical obstacles (d = 3) with Ds = D b . Lines: theoretical prediction, Eq.[START_REF] Burada | Diffusion in confined geometries[END_REF]. Squares: data obtained by numerical integration of Eqs. (12),[START_REF] Carusela | Computational model for membrane transporters. potential implications for cancer[END_REF]. Triangles: results of stochastic simulations of the stochastic process r(t) associated to the Fokker-Planck equations (1,2). Circles: Monte Carlo simulation data taken from[START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF], where the tracer particles are submitted to an exponentially decaying potential. The parameter δ was adjusted so that the stationary probability densities (9) in our description are the same as in the case of this potential (see SI for details).

FIG. 3 .

 3 FIG. 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/ √ 5) with u = 2πx/L. Symbols: results of the numerical integration of Eqs. (12) and (13) for kaa/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.

= 1 /

 1 (1 + (2δ/R -1)ϕ) and P st,0 s = δP st,0 b

5 r4must match F 7 b+

 57 D32) Next, the inner and outer solutions at fourth order vanish, f 4 b = f 4 s = F 4 b = 0. The inner solution f 5 b reads: f 5 b (r, θ) = [3 cos(θ) + 5 cos(3θ)] A D 0 [3 cos(θ) + 5 cos(3θ)] r3 matching F 0 b , f 5 s (θ) = [3 cos(θ) + 5 cos(3θ)] f 5 s . (D33) The constants A 5 and f 5 s are found by solving (D24) with i = 5 but their precise expressions are not necessary because π 0 dθ 2π 0 dφ cos(θ) sin(θ) [3 cos(θ) + 5 cos(3θ)] = 0 , (D34)

FIG. 4 .

 4 FIG. 4. Left:Comparison between Monte Carlo simulation data for dispersion in 3D arrays of obstacles (circles: data of Fig.9in the SI of Ref.[START_REF] Putzel | Nonmonotonic diffusion of particles among larger attractive crowding spheres[END_REF]), our theoretical prediction and the result (D45) of Ref.[START_REF] Edwards | Diffusion and convective dispersion through arrays of spheres with surface adsorption, diffusion, and unequal solute partitioning[END_REF]. Parameters: ϕ = 0.57, Ds = D b , γ = 0, and δ is identified using Eq. (D44) 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: ϕ = 0.2, Ds = D b , γ = 0 (blue) or γ = 0.1 (orange) and a time step dt = 10 -6 L 2 /D b .

1 h

 1 (D44)], and we note that a d = ds d . The free energy F as a function of x is thus given byF(x) = -1 β ln h(x) d-1 + (d -1)δh(x) d-2 + ln(s d-1 ) . (E33)Also, we can introduce the average diffusion constant at position x given byD a (x) = D s (d -1)δh(x) d-2 + D b h(x) d-1 (d -1)δh(x) d-2 + h(x) d-1 . (E34)The effective diffusion constant in one dimension in a periodic potential F(x) with periodic diffusion constant D a (x) is given by the Lifson-Jackson formula[START_REF] Lifson | On the self-diffusion of ions in a polyelectrolyte solution[END_REF] D e = 1 exp(-βF) exp(βF)D -1 a = d-1 + (d -1)δh d-2 (D b h d-1 + D s (d -1)δh d-2 )

= 2 ,

 2 k d fs -k a fb , ∂x = k d fs -k a fb + D b 1 + w 2 (r)a

  ) , r] = 0 , D b ∂ f 0 b ∂x [ w(r), r] = k d f 0 s (r) -k a f 0 b [ w(r), r] + D b . (F2)Solving these equations leads to:f 0 b (x, r) = A C (r) for |r| < 1, A P (r) + x for |r| > 1, (F3)where the functions A C and A P depend only on 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 (F2). 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 the general expression of D e yields D e = a→+∞ D b P st b,0 (V -V p ) , (F4)

Appendix E: Dispersion in slowly undulated channels (a L) -Proof of Eq. ( 18)

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 e s,x = e x + h (x)e y , so that

. Using the formulas of Section B, Eq. (12,