M Mangeat

T Guérin

D S Dean

Dispersion in two dimensional channels-the Fick-Jacobs approximation revisited

Keywords: Brownian motion, diusion, transport processes/heat transfer

We examine the dispersion of Brownian particles in a symmetric two dimensional channel, this classical problem has been widely studied in the literature using the so called Fick-Jacobs' approximation and its various improvements. Most studies rely on the reduction to an eective one dimensional diusion equation, here we derive an explicit formula for the diusion constant which avoids this reduction. Using this formula the eective diusion constant can be evaluated numerically without resorting to Brownian simulations. In addition, a perturbation theory can be developed in ε = h 0 /L where h 0 is the characteristic channel height and L the period. This perturbation theory confirms the results of Kalinay and Percus (2006 Phys. Rev. E 74 041203), based on the reduction, to one dimensional diusion are exact at least to O(ε 6 ).

Furthermore, we show how the Kalinay and Percus pseudo-linear approximation can be straightforwardly recovered. The approach proposed here can also be exploited to yield exact results in the limit ε → ∞, we show that here the diusion constant remains finite and show how the result can be obtained with a simple physical argument. Moreover, we show that the correction to the eective diusion constant is of order 1/ε and remarkably has some universal characteristics. Numerically we compare the analytic results obtained with exact numerical calculations for a number of interesting channel geometries.

Introduction

The dispersion properties of complex systems have a wide range of applications including solid state physics, soft matter physics, biology and chemical engineering [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF][START_REF] Gardiner | Handbook of Stochastic Methods for Physics[END_REF]. In many cases the eective transport properties, such as the late time diusion constant and eective drift are important quantities to study problems as varied as mixing [3][4][5], sorting [6], contaminant spreading [7] or reaction kinetics in complex systems [8]. The computation of such eective transport quantities is a longstanding problem in statistical physics and exploits techniques ranging from probability theory [9] through to quantum field theory [10].

Transport, notably diusion, in channel like geometries arises in a number of important physical contexts, ranging from biophysics through to the study of porous media [11][12][13][14]. Dispersion in periodic systems such as hard squares or disks placed on a lattice can also be mapped onto diusion in channel like systems [15]. Even if the underlying process in such problems is locally standard Brownian motion, the interplay between geometry and diusion gives rise to rich behavior, for instance the slowing down of diusive dispersion due to entropic trapping. In many contexts one is only interested in dispersion along the channel and the problem of computing the eective diusion constant parallel to the channel is one which has been intensively studied. The first results on diusion in channels are attributed to Jacobs in [START_REF] Jacobs | Diusion Processes[END_REF] where the first form of what is known as the Fick-Jacobs' (FJ) approximation was derived. Subsequent works [START_REF] Zwanzig | [END_REF][18][19][20][21][22][23][24] involved improvements to the basic FJ approximation for both two dimensional channels and axisymmetric channels in three dimensions. These studies are based on the eective reduction of the problem of diusion in the channel, to a one dimensional problem where a particle diuses in an entropically generated potential, simply related to the local channel width or area, but with an eective local diusion constant. Once this reduction is carried out the eective diusion constant can be computed using exact results for the diusion constant of one dimensional systems [25][26][27].

One of the most successful approaches to the problem of diusion in a channel was proposed by Kalinay and Percus in [19,20]. Here they express the eective local diusion constant as an expansion in powers of a parameter ε proportional to the ratio of the typical relaxation times in the perpendicular and parallel direction. Later studies [28] checked that such expansions lead to expressions of the global eective diusivity that are exact at least to order O(ε 4 ). From a probabilistic point of view this latter result is somewhat surprising as the eective process parallel to the channel can clearly never be Markovian and described by an eective one dimensional diusion equation. Mathematically similar problems, where a higher dimensional Markovian stochastic processes is reduced to a lower dimensional one, occur throughout the physics literature, notably when partially-damped Langevin equations (describing the velocity and position of a particle) are reduced to over-damped Langevin equations (describing only the position). A similar problem arises for a particle carrying a dipole in a spatially varying electric field, here the Langevin equations describing the dipole moment and the position can be reduced to a Langevin equation for just the position. In the class of problems mentioned above, a mathematical process of adiabatic elimination, which amounts to a perturbation expansion in the time scale for the relaxation of the quick variable (velocity or dipole moment for instance) divided by that of the slow position variable can be applied. In most cases the eective process remains Markovian only to first order in perturbation theory [START_REF] Risken | The Fokker-Planck Equation 3rd edn[END_REF][START_REF] Sancho | [END_REF][31]. In this paper we neglect hydrodynamic eects, these can become important for fluid filled channels [32] when the particle size becomes comparable to the channel size. This is due to the spatial variation of the local diusion tensor as one approaches the walls channel. Here we restrict our attention to the dispersion of a point-like particle and we can thus ignore these hydrodynamic eects.

In this paper, we analyze in section 2 the FJ problem starting from an exact formula for the longitudinal diusion constant. The formula is derived here using a physical argument for channels in any dimension, and agrees with general Kubo type formulas derived for advection diusion equations in periodic systems [34,35] and with the macroscopic transport theory developed by Brenner [START_REF] Brenner | Macrotransport Theory[END_REF]. Then, for two-dimensional channels, we present in section 3 a very compact analysis which allows us to show in section 4 that the results of [20] are exact to O(ε 6 ), whereas to date they had been checked only up to O(ε 4 ) by using Brenner's theory [28]. Given what we have stated about adiabatic approximations above, it is very interesting that the method of [20] stays correct up to such a high order. We obtain this result within a formalism that uses analytic functions in a complex plane, and we show in section 5 how the pseudo-linear approximation of Kalinay and Percus can be recovered. This approximation is based on neglecting terms W (x) (where W is the local channel width) and higher order derivatives. In section 6, we show how the diusion constant can be computed exactly for a class of two dimensional channels which can be appropriately parameterized in the complex plane, this reproduces in a very rapid manner the classical result of Zwanzig [START_REF] Zwanzig | [END_REF]. We also show in section 7 how the same formalism can be applied to the limit ε → ∞, the exact opposite of the cases usually studied in the literature and in which all approaches based on a dimensional reduction fail. In this limit ε → ∞, we find the asymptotic value of the diusion constant (which is non-zero) and we also compute the leading order correction, which is O(1/ε), and which has a prefactor that has a remarkably universal character. A short description of this result was given in [38], the derivation presented here is slightly dierent, as we concentrate here on two dimensional channels, and relies more heavily on complex analysis. Finally, in section 8 we show how the formalism based on the Kubo formula for the diusion constant can be used to numerically study a wide range of channel geometries, thus avoiding the statistical errors and longer times associated with Brownian simulations. We compare some of the standard variants of the FJ approximation with these results for a number of channel geometries, as well as the high order perturbation theory and large ε results derived here.

General equations for diusion in channels

Model

In this section, we consider diusion in a portion of space located between two surfaces of equations y = h 1 (x) and y = h 2 (x). Here x represents the coordinates in the hyperplane parallel to the channel, and y is the coordinate in the perpendicular direction. We consider the diusion of a tracer with isotropic diusion in the channel, with diusion constant D, and subject to a uniform drift velocity parallel to the channel u || . The probability density function (pdf) for the tracer thus obeys the advection-diusion equation

∂p ∂t (x, y; t) = ∇ • D∇p(x, y; t) -u || p(x, y; t) . ( 1
)
with no-flux boundary condition at the interface. We now define the marginal pdf for the coordinate x as

p * (x; t) = h 2 (x) h 1 (x)
dy p(x, y; t).

(

) 2
The basis of FJ type approximations is to find an eective diusion equation for the marginal pdf in the case where it describes a one dimensional process. In what follows, we will (i) derive an exact equation for p * (x, t) to make the hypotheses of the FJ approximation more explicit, and (ii) we will derive an exact equation for the eective diusivity.

An exact equation for the evolution of the marginal distribution

Here, we shall derive an exact equation for p * , which makes more explicit the approximations of the FJ approach. To proceed we take the time derivative of equation (2) and use the diusion equation (1) to obtain (3) The no-flux boundary condition at each surface is given explicitly as at the surface y = h α (x) (for α = 1, 2):

∂p * ∂t (x; t) = h1(x) h2(x)
D ∂p ∂y (x, h α (x); t) = [D∇ || p(x, y; t) -u || p(x, y; t)]| y=hα(x) • ∇ || h α (x), (4)
where we have used the decomposition of the gradient operator as (5) Now using the fact that u || is uniform we can write (6) and we also find

∇ || • [u || p * (x; t)] = ∇ || • h2(x) h1(x) dy u || p(x, y; t) = ∇ || h 2 (x) • u || p(x, h 2 (x); t) -∇ || h 1 (x) • u || p(x, h 1 (x); t) + h2(x) h1(x) dy∇ || • [u || p(x, y; t)],
∇ 2 || p * (x; t) = ∇ 2 || h2(x) h1(x) dy p(x, y; t) = ∇ || • [ p(x, h 2 (x); t)∇ || h 2 (x)] + ∇ || h 2 (x) • ∇ || p(x, y; t)| y=h2(x) -∇ || • [ p(x, h 1 (x); t)∇ || h 1 (x)] -∇ || h 1 (x) • ∇ || p(x, y; t))| y=h1(x) + h(x) h2(x)
dy ∇ 2 || p(x, y; t).

(7) Using these results equation (5) simplifies remarkably to give

∂p * ∂t (x; t) = D∇ 2 || p * (x; t) -∇ || • [u || p * (x; t)] -D∇ || • p(x, h 2 (x); t)∇ || h 2 (x) -p(x, h 1 (x); t)∇ || h 1 (x) .
(8) The above equation is exact and is the main result of this section. Even if it is not a closed equation for p * , we see that the evolution of the marginal pdf p * involves only the value of the full probability density at the surfaces y = h α (x), not in the whole channel.

The basic Fick-Jacobs' approximation

The basic FJ approximation is obtained by assuming that p(x, y, t) p 0 (x, t), i.e. assuming that it is independent of y, and thus one is assuming that the pdf equilibrates in the coordinate y. This means that p * (x; t)

h 2 (x) h 1 (x)
dy p 0 (x; t) = H(x) p 0 (x; t), (9) where

H(x) = h 2 (x) -h 1 (x)
is the total channel width at the in-plane coordinate x.

The basic FJ approximation then gives p(x, h 1 (x); t) p(x, h 2 (x); t) p * (x; t) H(x) (10) and so the eective diusion equation for p * (x; t) becomes

∂p * ∂t (x; t) = D∇ 2 || p * (x; t) -∇ || • u || p * (x; t) -D∇ || • p * (x; t)∇ || ln[H(x)] . ( 11
)
The physical validity of this approximation is argued to hold when the process relaxes in the y direction very quickly with respect to its motion in the direction parallel to the channel; it should hold when H(x) L where L is the length scale over which the total channel width varies. An alternative limit where it holds is when one takes the diusion to be anisotropic, having a (isotropic) diusion constant D || in the plane and a dierent diusion constant D ⊥ in the perpendicular direction. The FJ approximation then holds in the limit ε = D || /D ⊥ 1, and the problem can be analyzed perturbatively using ε as the small parameter. Of course the approximation that p(x, y; t) is independent of y clearly does not satisfy the no flux boundary condition in equation (4) when ∇ || h α = 0. An important success of the FJ approximation is that it gives the correct equilibrium distribution in the steady state. On a unit cell Ω of the channel (which is repeated to form the infinite channel) the equilibrium distribution is clearly uniform within the channel and thus the marginal pdf for the in plane coordinate x is

P s (x) = H(x) Ω dx H(x ) . ( 12
)
This means that the x in-plane coordinate of the tracer particle is subjected to an eective entropic potential

φ ent (x) = -k B T ln(H(x)). ( 13
)
For eectively two dimensional systems (that is to say two dimensional channels or three dimensional ones with rotational symmetry) a number of authors have suggested improvements to the Fick-Jacobs approximation in the form an eective in plane diusion constant D(x) which depends on the local height functions (in general both h 2 and h 1 and not just their dierence). Clearly if an underlying one dimensional diusion equation were to exist it would have to have the steady state distribution given by equation (12) and thus have the general form, in the absence of external drift,

∂p * ∂t (x; t) = ∂ ∂x D(x) ∂p * ∂x (x; t) + β ∂φ ent ∂x (x) p * (x; t) , ( 14
)
with β =1/k B T . The eective diusion constant for the eective one dimensional prob-lem can then be extracted using the exact formula for one dimensional diusion [25]

D e = 1 D -1 (x) exp (βφ ent (x)) s exp (-βφ ent (x)) s , ( 15
)
where the angled bracket indicate the spatial average over one period of the system, i.e.

f s = 1 L L 0 dx f(x). ( 16
)
These approaches all agree in an expansion based on ε = D ⊥ /D || being small but dier at higher orders. However, as mentioned in the introduction, this approach clearly cannot be valid exactly as the process projected onto the in-plane coordinates is evidently not Markovian and thus cannot be described by a diusion equation with local transport coecients. As mentioned in the introduction, a similar phenomenon appears in the theory of stochastic processes where a fast variable is eliminated from a stochastic dierential equation using an adiabatic approximation scheme [START_REF] Risken | The Fokker-Planck Equation 3rd edn[END_REF][START_REF] Sancho | [END_REF][31]. It is generally found that the process remains Markovian at the first order of correction but does not remain so when higher order corrections are included.

General formula for the diusion constant

In this paper we will use an exact formula for the diusion constant, of the process in the plane of the channel, which does not use a reduction of the problem to a lower dimensional one. The formulas for the diusion constant can be directly written down using a general theory for diusion in periodic systems developed in [34,35]. These formulas are formally exact-however, they require the solution of partial dierential equations which are not analytically soluble in two dimensions and higher. Nevertheless, the relevant equations can be treated using approximation methods for narrow channel geometries which are similar to those which are employed to derive improvements to the FJ approximation and which are used as the basis for lubrication theory in fluid dynamics [33].

From now on, we consider that h 1 and h 2 are spatially periodic, and of the same periodicity. Rather than simply quoting the results of [34,35], we provide a quick re derivation of the result, ideally suited to the channel geometry and based on simple physical arguments. To begin with we compete the average drift of a tracer in the in plane coordinates,

d X i dt = dxdy x i ∂p ∂t (x, y; t) = dx x i ∂p * ∂t (x, t), ( 17
)
and we now use equation (8) and integration by parts to give an eective drift

d X dt = u || + D dx p (x, h 2 (x); t) ∇ || h 2 (x) -p (x, h 1 (x); t) ∇ || h 1 (x) . ( 18
)
In the steady state we can replace p (x, h 2 (x); t) by their steady state values on the unit cell Ω, denoted by P s (x, y) and use periodicity to restrict the x integral over the in-plane coordinates of Ω denoted by Ω || , to find

d X dt = u || + D Ω || dx P s (x, h 2 (x)) ∇ || h 2 (x) -P s (x, h 1 (x)) ∇ || h 1 (x) . ( 19
)
Now we suppose that the drift u || is induced by an infinitesimal applied force F || , i.e. u || = DβF || . This local constant applied force will give rise to an eective late time, or asymptotic, drift which is defined by

V ||i = d X i /dt.
The Stokes Einstein relation then gives

V ||i = βD ij e F || j = D ij e D u || j , ( 20
)
where the Einstein convention on index summation is used, and D ij e is the eective diusion tensor in the in-plane coordinates, given for u || = 0 by

D ij e = lim t→∞ x i (t)x j (t) 2t . ( 21
)
The steady state distribution obeys

∇ • D∇P s (x, y) -u || P s (x, y) = 0, (22)
and then assuming that u ⊥ = 0 and u || is small we can write

P s (x, y) = 1 Ω + u || j ∂P s ∂u || j (x, y), ( 23
)
where we have also used Ω to denote the volume of the unit cell Ω. We write

f i = D ∂Ps ∂u ||i
and so the Stokes Einstein formula gives

D ij e D = δ ij + Ω || dx f j (x, h 2 (x)) ∇ ||i h 2 (x) -f j (x, h 1 (x)) ∇ ||i h 1 (x) . ( 24
)
Dierentiating the steady state equation (22) with respect to u ||i and then setting u ||i = 0 yields the equation obeyed by the vector field f i :

∇ 2 f i = 0. ( 25
)
Furthermore, f i (x, y) must be periodic in the in plane coordinates x and the normalization of the steady state distribution over Ω implies that

Ω dydx f i (x, y) = 0, (26)
while the no flux boundary conditions at the channel surface imposes

n α • -∇f i (x, y) | y=hα(x) + e i Ω = 0, (27)
where e i denotes the unit vector in the direction i in the plane parallel to the channel. The above condition can be rewritten as

-∇ || h α (x) • ∇ || f i (x, y) | y=hα(x) + ∂f i ∂y (x, h α (x)) = - ∇ ||i h α (x) Ω . ( 28
)
To summarize, an exact formula for the eective diusivity is given by equation (24), where the auxiliary function f i is computed by solving Laplace's equation (25) with boundary conditions (28). The above results are in full accordance with the general theory of [34,35] and with the Brenner macro transport theory [START_REF] Brenner | Macrotransport Theory[END_REF].

Formulation of the problem with analytic functions for two dimensional symmetric channels

Most studies in the literature have considered the case where the channel is two dimensional. Here we will denote the coordinate along the channel as x and, as earlier, use y to denote the height coordinate. We will consider channels which are symmetric about the axis y = 0, that is to say we write h(x) = h 2 (x) = -h 1 (x). We define L as the period of the channel, so that h(x

+ L) = h(x).
Here there is only one variable f 1 = f corresponding to the x direction it obeys Laplace's equation

∇ 2 f (x, y) = 0. ( 29
)
The solution of Laplace's equation can be written using the theory of analytic functions as the real or imaginary part (or a linear combination of both) of

F (x, y) = W 1 (z) + W 2 (z), (30)
where z = x + iy and z = x -iy denotes its complex conjugate. For symmetric channels, the symmetry about y = 0 implies that F (x, y) = F (x, -y) and thus one can easily see that

W 1 (z) = W 2 (z) + c, ( 31
)
where c is a constant. However, c can simply be incorporated into the definition of W 1 so we can use the form

F (x, y) = W 1 (z) + W 1 (z). (32)
Since f (x, y) is real, we can assume without loss of generality that the coecients of the Laurent expansion of W 1 are real, so that W 1 (z) = W 1 (z). Making a convenient rescaling in terms of the unit channel area Ω we thus see that, for a symmetric channel, one can look for solutions of the form

f (x, y) = 1 2Ω [w(x + iy) + w(x -iy)], ( 33
)
where w is analytic such that w(z) = w(z) and is periodic in x (with period L). If we write w in terms of its real and imaginary parts, w(x + iy) = u(x, y) + iv(x, y), the symmetry y → -y implies that

u(x, -y) = u(x, y). ( 34
)
This means, looking at an expansion in terms of y at fixed x, that ∂u/∂y is an odd function in y. However, the Cauchy-Riemann equation ∂u ∂y (x, y) = -∂v ∂x (x, y), (35) shows that ∂v/∂x is odd in y. However, at y = 0 this equation implies, as ∂u ∂y | y=0 = 0, that v(x, 0) is a constant (v 0 ). Consequently, we have, integrating equation (35)

over x, that v(x, y) = v 0 + V (x, y) ( 36
)
where V (x, -y) = -V (x, y). However there is nothing in the problem that allows the determination of v 0 , furthermore one can show that it does not aect the value of the diusion constant, and we thus chose v 0 = 0 and so v(x, -y) = -v(x, y).

Inserting the ansatz in equation (33) into the boundary condition (28) provides an additional condition for the value of w at the channel boundary:

ih (x)[w (x + ih(x)) + w (x -ih(x))] + [w (x + ih(x)) -w (x -ih(x))] = i2h (x),
(37) and this can be integrated with respect to x to yield

w (x + ih(x)) -w (x -ih(x)) = 2ih(x) -i2C, ( 38
)
where C is a constant to be determined. Finally, the eective diusivity is given by

D e D = 1 + 2 L 0 dxf (x, h(x))h (x), (39)
and so within this formalism, the diusion constant is simply given by

D e D = 1 + 1 Ω L 0 dx[w (x + ih(x)) + w (x -ih(x))]h (x). ( 40
)
However, we show in appendix that the above equation can be simplified, leading to

D e D = C h s , ( 41
)
where

h s = 1 L L 0 dx h(x), ( 42
)
is the channel height spatially averaged over one period. In the following section we present an alternative proof based on a perturbation expansion.

We have thus reformulated the diusivity problem for two dimensional symmetric channels into a problem of finding an analytic function w(z) in the complex plane z = x + iy, with the condition (38) at the channel boundary. Remarkably, once the constant C is determined in (38), there is no more integration to be done to determine D e which is directly proportional to this constant (see equation (41)). The use of analytic functions to analyze FJ type problems in two dimensional channels has of course been exploited in a number of studies [39][40][41]; however, the present study shows that this method can be exploited much further.

Perturbation expansion for the eective diusivity

Recovering the basic Fick-Jacobs' approximation

We now show how our complex formalism can be used to recover the basic FJ approximation for the eective diusivity. To proceed we assume that the length of the channel is much larger than the height and thus expand equation (38) for small h which gives

w (x)h(x) = h(x) -C. ( 43
)
Integrating and using the periodicity of w we find

w(x) = x -h -1 -1 s x 0 dx h(x ) + C ( 44
)
where

h -1 s = 1 L L 0 dx h(x ) , ( 45
)
is the spatial average of h -1 over the unit cell. The term C is determined from the normalization condition but we do not need to evaluate it gives zero contribution to the eective diusivity as can be easily seen in equation (40). Keeping the leading order in h in equation (40) now gives

D e D 1 - 2 Ω L 0 dx w (x)h(x). ( 46
)
Now, noting that

Ω = 2 L 0 dx h(x) we find D e D ≈ D FJ D = 1 h -1 s h s , ( 47
)
this is exactly the same result as given by the basic FJ formula (that is to say with no modifications to the local diusion constant).

Perturbation theory in channel width

We now consider a channel height of the form

h(x) = εζ(x), ( 48
)
which we will analyze here in the limit of small ε. We now show that our complex formalism can be systematically improved to obtain an expansion of the eective diusivity in powers of the small parameter ε in a very compact way.

The whole perturbation analysis for the function w(x) is based on equation (38) which can be written as

∞ n=0 (-1) n ε 2n+1 w (2n+1) (x)ζ 2n+1 (x) (2n + 1)! = εζ(x) -C, ( 49
)
where we use the notation w (n) = d n w/dx n . The above equation is compatible with an expansion of w in even powers of ε of the form

w(x) = ∞ m=0 ε 2m w 2m (x), (50)
and the constant C of the form

C = ε ∞ k=0 ε 2k C 2k . ( 51
)
The boundary condition equation (38), expressed as a perturbation series in ε, now becomes for k > 0,

∞ n,m=0 (-1) n ε 2n+2m w (2n+1) 2m (x)ζ 2n+1 (x) (2n + 1)! = ζ(x) - ∞ k=0 ε 2k C 2k . ( 52
)
We can compute directly using equation (39), the expansion in ε of the eective diusivity, finding

D e D = 1 - 1 ζ s L 0 dx ∞ n=0 (-1) n ε 2n w (2n+1) (x)ζ 2n+1 (x) (2n + 1)! , ( 53
)
after an integration by parts. It follows from the relation equation ( 49) that the constants C i directly provide the terms of the perturbative expansion of the eective diusion constant, which reads

D e D = C ε ζ s = 1 ζ s ∞ k=0 C 2k ε 2k , ( 54
)
thus providing an alternative derivation of equation (41) to that derived in appendix.

Let us write explicitly the first orders in ε k of equation (52). The leading order term in ε gives the basic FJ approximation found above

ζ(x)w (1) 0 (x) = ζ(x) -C 0 , ( 55
)
while the higher order corrections give for k > 0,

k n=0 (-1) n w (2n+1) 2(k-n) (x)ζ 2n+1 (x) (2n + 1)! = -C 2k . ( 56
)
The next three orders in perturbation theory provide the equations satisfied by w 2 , w 4 , w 6 and read:

- 1 3! w (3) 0 (x)ζ 3 (x) + w (1) 2 (x)ζ(x) = -C 2 , ( 57
)
for k = 1, 1 5! w (5) 0 (x)ζ 5 (x) - 1 3! w (3) 2 (x)ζ 3 (x) + w (1) 4 (x)ζ(x) = -C 4 , ( 58
)
for k = 2 and - 1 7! w (7) 0 (x)ζ(x) 7 + 1 5! w (5) 2 (x)ζ(x) 5 - 1 3! w (3) 4 (x)ζ(x) 3 + w (1) 6 (x)ζ(x) = -C 6 , ( 59
)
for k = 3.
We now show how to compute the constants C 2k . For k = 0, the value of C 0 is found by imposing periodic boundary conditions in equation (55), which gives

C 0 = ζ -1 -1 s . ( 60
)
Inserting into equation (54) one sees that one of course recovers the basic form of the FJ approximation, the solution for w 0 is in this notation

w (1) 0 (x) = 1 - 1 ζ -1 s ζ(x) . ( 61
)
Substituting this solution into equation (57) now gives

w (1) 2 (x) = - C 2 ζ(x) + ζ -1 -1 s 6 ζ (2) (x) -2 ζ (1) (x) 2 ζ(x) ( 62
)
and the periodicity condition w 2 (0) = w 2 (1) now gives

C 2 = - 1 3
ζ 2 /ζ s ζ -1 2 s . ( 63
)
Substituting the solution for w

2 into equation (58) yields

w (1) 4 (x) = - C 4 ζ(x) + 1 6 ζ 2 (x) - C 2 ζ(x) + C 0 6 ζ (2) (x) -2 ζ (1) (x) 2 ζ(x) (2)
+ C 0 120 ζ 4 (x) 1 ζ(x) (4)
.

(64) The periodicity of w 4 then gives, after some integration by parts,

C 4 C 0 = C 2 2 C 2 0 + C 0 45 4 ζ 4 /ζ s + ζ 2 ζ s , ( 65
)
which can be rearranged as

C 4 = ζ 2 /ζ 2 s 9 ζ -1 3 s + 1 45
4 ζ 4 /ζ s + ζ 2 ζ s ζ -1 2 s . ( 66
)
Using the above expressions for C 0 , C 2 et C 4 we find to O(ε 4 ):

D e D = 1 ζ s ζ -1 s 1 -ε 2 ζ 2 /ζ s 3 ζ -1 s + ε 4 ζ 2 /ζ 2 s 9 ζ -1 2 s + 4 ζ 4 /ζ s 45 ζ -1 s + ζ 2 ζ s 45 ζ -1 s . ( 67
)
This agrees with the results of [20,28]. In [20] the result was obtained by a formalism which derives an eective one dimensional diusion constant. However, as pointed out before, the reduction to an eective one dimensional problem is not obvious as the resulting one dimensional process is non-Markovian. In [28] the eective macroscopic transport theory of Brenner [START_REF] Brenner | Macrotransport Theory[END_REF] was used. The expression for the diusion constant in this theory can actually be shown to agree with the Kubo formula derived here. An interesting point made in [28] is that the result of [20] is recovered at orders up to and including ε 4 .

Here we exploit the compact nature of the analysis here to derive the term of O(ε 6 ) and show that it still agrees with [20]. Substituting the solution for w into equation (59) and imposing the periodicity of w 6 gives the value of C 6 as

C 6 C 0 = - C 4 3 ζ 2 /ζ s + C 2 45 4ζ 4 /ζ + ζ 2 ζ s - C 0 945 44ζ 6 /ζ + 5ζ 2 ζ 3 -60ζ 4 ζ + 9ζζ 2 ζ 2 -4ζζ 3 ζ + 2ζ 2 ζ ζ ζ -2ζ 3 ζ ζ (4) + 4ζ 2 ζ 2 ζ (4) s .
(68) Using integration by parts to minimize the highest derivative of ζ appearing in the final expressions we obtain

C 6 C 0 = C 3 2 C 3 0 + 2C 4 C 2 C 2 0 - C 0 945 44 ζ 6 /ζ s + 5 ζ 2 ζ 3 s + 45 ζ 3 ζ 2 s + 2 ζ 3 ζ 2 s , ( 69
)
which can be rearranged as

C 6 = - 1 27
ζ 2 /ζ 3 s ζ -1 4 s - 2 135
ζ 2 /ζ s 4ζ 4 /ζ + ζ 2 ζ s ζ -1 3 s - 1 945
44 ζ 6 /ζ s + 5 ζ 2 ζ 3 s + 45 ζ 3 ζ 2 s + 2 ζ 3 ζ 2 s ζ -1 2 s (70)
The final expression for the diusion constant to O(ε 6 ) can also be conveniently written as

D e D = ζ -1 s ζ -1 s + ε 2 3 ζ 2 /ζ s - ε 4 45 4 ζ 4 /ζ s + ζ 2 ζ s + ε 6 945 44 ζ 6 /ζ s + 5 ζ 2 ζ 3 s + 45 ζ 3 ζ 2 s + 2 ζ 3 ζ 2 s -1 , ( 71
)
which is also true at order O(ε 6 ) and after some work one can show that this agrees with the result of [20] at O(ε 6 ).

The pseudo-linear channel approximation

In [20] the perturbation theory was partially resumed by assuming that derivatives ζ and of higher order vanish. Here we show how this re-summation can be derived within the formalism developed here in a very compact manner. We start from the equation (38) and make the ansatz

w(x + iy) = x + iy - 1 φ s x+iy 0 dξφ(ξ), (72)
where

φ s = 1 0 φ(x)dx (73)
is the integral of φ along the line y = 0. In addition, we assume that φ is analytic and so the integral in the second term on the right hand side is independent of the integration path between 0 and x + iy. We also assume that φ is of period 1 in the variable x, i.e. φ(x + 1 + iy) = φ(x + iy), so we are working in units where the period L = 1 to simplify notation.

We now note that the function w so constructed is analytic and is periodic in x with period 1. To see this latter point we write

w(x + 1 + iy) = x + 1 + iy - 1 φ s x+1+iy 0 dξφ(ξ). ( 74
)
The contour integral on the right hand side can be decomposed into an integral along the real axis between 0 and x + 1 followed by an integral parallel to the imaginary axis between x + 1 and x + 1 + iy so we can write

w(x + 1 + iy) = x + 1 + iy - 1 φ s x+1 0 dx φ(x ) - i φ s y 0 dy φ(x + 1 + iy ) (75)
However, by the periodicity of φ we can write

1 φ s x+1 0 dx φ(x ) = 1 + 1 φ s x 0 dx φ(x ) (76)
which once inserted into (74) proves the desired periodicity, i.e. w(x + 1 + iy) = w(x + iy).

We now substitute this ansatz into the boundary condition (38) which becomes

1 φ s x+ih(x) 0 dξφ(ξ) - x-ih(x) 0 dξφ(ξ) = 2iC. ( 78
)
We now note that the overall contour integral on the left hand side goes over the line between x -ih(x) and 0 then the line from 0 to x + ih(x). Using the analytic property of φ we deform the contour to the straight line between x -ih(x) and x + ih(x) which is parallel to the imaginary axis so we can write

1 φ s 1 -1 ds h(x) φ (x + ish(x)) = 2C. ( 79
)
It is easy to verify that if we write h(x) = εζ(x) the above generates the same perturbation theory as was found in the previous section. However, the starting point of the perturbation theory can be changed by writing

φ(ξ) = ψ(ξ) h(ξ) , ( 80
)
which leads to

1 ψ h s 1 -1 ds ζ(x) ψ (x + isεζ(x)) ζ (x + isεζ(x)) = 2C. ( 81
)
At lowest order in ε it is easy to verify using equation (54) that this recovers the basic FJ approximation as the solution is ψ(x) = ψ 0 where ψ 0 is constant.

To recover the pseudo-linear approximation of [20], we proceed as follows. First, we ignore the terms ζ (x) and higher orders in Taylor expanding the denominator of the integrand on the right hand side of equation (81), so that we approximate

ζ(x + isεζ(x)) ζ(x) + ζ (x) isεζ(x).
(82)

Next, we also write

ψ (x + isεζ(x)) ≈ ψ(x), ( 83
)
later we will verify that this is consistent with neglecting terms h and higher derivatives. With these approximations, equation (81) becomes

ψ(x) ψ h s 1 -1 ds 1 + isεζ (x) = 2C. ( 84
)
Performing the integral in the above expression leads to

ψ(x) = C ψ h s ε|ζ (x)| tan -1 (ε|ζ (x)|) . ( 85
)
where tan -1 is the arctangent function. Here we see that ψ (x) ∼ h (x) and thus the approximation made in (83) is consistent. Finally, multiplying both sides of equation (85) by 1/h = 1/(εζ) and integrating over [0, 1] we find

C = |h (x)| h(x) tan -1 (|h (x)|) -1 s , ( 86
)
and thus equation (41) gives

D e D = 1
h s h (x) h(x) tan -1 (h (x)) s . ( 87
)
Within their systematic reduction to a one dimensional diusion equation the authors of [20] show that the eective one dimensional diusion constant that should be used in equation (15) is

D(x) = D tan -1 (h (x)) h (x) , ( 88
)
and hence by using the Lifson-Jackson formula (15) we recover exactly the same formula (87) under the same approximation.

Generating exact results

Here we re derive in a very direct manner a classical exact result due to Zwanzig [START_REF] Zwanzig | [END_REF] for channels which can be appropriately parameterized in another coordinate system. Consider the transformation

z = G(w) ( 89
)
where z = x + iy and w = u + iv where G is an analytic function. Importantly, the upper and lower channel boundaries are defined to be given by v = ±V where V is constant. Of course the function G must be chosen to be periodic in u and should be single valued. Equation (78) is trivially written as

1 φ s x+ih(x) x-ih(x) dξ φ(ξ) = 2iC, ( 90
)
where the integral is a contour integral over any finite contour joining the lower to upper limits indicated on the integral. Now in terms of the variables w this becomes

1 φ s u+iV u-iV dw dG dw φ(w) = 2iC. ( 91
)
The solution to this equation (up to a multiplicative constant) is now clearly given by

φ(w) = dG dw -1 , ( 92
)
which gives

C = V dG dw -1 s , ( 93
)
and equation (41) now becomes

D e D = V dG dw -1 s h s . ( 94
)
This can be rewritten by noticing that h s = Ω/2 is half the area of the unit cell and so we can write in the coordinate system given by w

h s = 1 2 1 0 du V -V dv dG dw 2 . ( 95
)
For a symmetric channel we can make the line y = 0 coincide with the line v = 0, so G is real for v = 0, dG dw

-1 s = 1 0 dx dG dw v=0 -1 = 1 0 duG (u)[G (u)] -1 = 1, ( 96
)
where we have used the notation dG dw | v=0 = G (u). Putting this altogether then gives

D e D = 2 V 1 0 du V -V dv| dG dw | 2 . ( 97
)
Up to dierences in notation the above is identical to the result of Zwanzig [START_REF] Zwanzig | [END_REF].

The large ε limit

Here we can consider the case of very wide channels where we again write h(x) = εζ(x) but we assume that ε → ∞. Before analyzing this problem via the formalism developed above we propose a simple physical argument which gives a bound on the diusion constant. Consider a symmetric channel whose height function about the center of the channel achieves its minimum value at some point x min where it takes the value h min .

Without loss of generality we choose x min = -1/2 ≡ 1/2. Now while the particle is in the region Ω = [-1/2, 1/2] × [-h min , h min ] the particle diuses without hindrance from the sides of the channel, so we can write

d dt X 2 = 2D, (98)
however, outside of this region the diusion is slowed down. If we denote by t Ω the total time spent in Ω we can thus write

X 2 2Dt Ω , ( 99
)
because time spent in the region where the walls can interfere with the diusion will nonetheless contribute to diusion-for example consider small excursions into this region which come back into Ω having moved along the channel.

As we know the equilibrium distribution in the channel modulo the period is uniform, ergodicity then implies that

t Ω t = |Ω | |Ω| = h min h s (100)
which gives the bound

X 2 2D h min h s t ( 101
)
and thus we have the bound on the eective diusion constant along the channel

D e D h min h s = ζ min ζ s . ( 102
)
It is also clear that D e D and thus we have the two sided bound

1 D e D ζ min ζ s , ( 103
)
and this simple result shows one that the system has a finite diusion coecient even in the limit ε → ∞. In [42] diusion on comb-like locally one dimensional structures was studied. The comb consists of a one dimensional backbone along which the diusion is measured with teeth perpendicular to the backbone, excursions into the teeth slow down diusion along the background. In this case the lower bound in equation (103) becomes exact as no diusion can occur in the teeth due to their one dimensional character and given they are perpendicular to the backbone. Other problems of biased diusion in tubes with dead-ends are studied in [43][44][45].

To go beyond this analysis we consider a channel where the points x = ±1/2 correspond to the minimum of the height of the channel. The height at x = ±1/2 is, as above, denoted by h min . To start with we define a function w(z) via

w(z) = z + w(z), (104)
where w(z) is the solution sought in equation (33). The boundary equation ( 38) contain now becomes

w (x + ih(x)) -w (x -ih(x)) = -i2C, ( 105
)
however, the function w is no longer periodic and we find

w(x + 1 + iy) -w(x + iy) = -1. ( 106
)
Now we make the transformation ξ = exp [-iπ(z -ih min )]. This conformal mapping is made explicit in the simple case where Ω is the septate channel shown in figure 1 and where we denote by h max the maximal channel height. The domain shown in figure 1(a) also represents the region in which the original function w(z) is analytic. In figure 1(b) we show the same region but include the branch cuts over which the function w(z) jumps. In this case the domain Ω is mapped into that shown on the figure 1(c), where the radius of the outer half-circle C γ is given by γ = exp [π(h max -h min )] and that of the inner circle

C λ is λ = exp [-π(h max + h min )].
The crucial point is that in the limit ε → ∞ the radius γ becomes infinite while the radius λ goes to zero (we recall that h is proportional to ε). For a general channel the points on the upper surface are mapped to a point at infinity while those on the lower surface are mapped to 0. We therefore argue that the obtained results will be valid for large ε for arbitrary channel shapes, since the details of the shape profile are sent to infinity after the mapping to the ξ space.

The resulting Laplace problem on the new domain corresponds to a classic two dimensional electrostatic problem, one has two charged plates, where the potential is φ = 1/2 on an upper plate, between (0, 0) and (0, 1), and φ = -1/2 on the lower plate, between (0, 0) and (0, -1). The potential vanishes at infinity. The solution is given by w(ξ) = i π arcsinh(1/ξ). The solution w(z) found here can include an arbitrary, undetermined, imaginary constant, corresponding to the unimportant constant v 0 in equation [START_REF] Brenner | Macrotransport Theory[END_REF]. Taking the logarithmic representation of arcsinh, we find

w(x + iy) = i π ln 1 + 1 + exp (-2πix + 2π(y -h min )) + i C. ( 107
)
where C is a constant. We note that, for y h min and h min 1, we have w(x) = i ln(2)/π + i C, and we recall that the choice v 0 = 0 means that the imaginary part vanishes on the central line y = 0, this fixes the value of C = -ln(2)/π. The above equation is valid for y > 0. For y < 0 the solution can be determined by symmetry to be w(z) = w(z).

Next we note that w(x + iy) -w(x -iy) in equation ( 38) is simply 2iIm[w(x + iy)], so that we can rewrite this equation and identify the constant C, considering the limit h(x)

1,

C = h min + ln 2 π ( 108
)
and finally we use equation (41) to find

D e D = ζ min ζ s + 1 ε ln 2 π ζ s . ( 109
)
and we thus recover within our complex formalism the result of [38]. This result clearly agrees with the ergodic argument given earlier when ε → ∞. We see that for finite ε the diusion is faster for the reasons evoked above. However, this result is rather intriguing, firstly the fact that the correction scales as 1/ε is far from obvious. Another remarkable feature is that this correction depends only on the average height function through 1/ ζ s and moreover the associated prefactor ln(2)/π is universal. These analytical results are valid for any channel shape in the limit of wide channels, as will be confirmed in a number of numerical tests in the following section.

Numerical results

In this section we will show how the diusion constant for periodic channels may be computed numerically. We use a standard finite element partial dierential equation solver to (i) solve Laplace's equation equation (25) with the integral constraint equation (26) and the boundary condition in equation ( 27) and (ii) carryout the subsequent evaluation of the eective diusion constant via equation (24). The results obtained are used to test some of the analytical results derived here and elsewhere in the literature. A jump of +1 in the solution occurs between the two curved lines.

The exact result of Zwanzig

In [START_REF] Zwanzig | [END_REF] Zwanzig considered a symmetric channel parameterized using the curvilinear coordinate system defined by

z = w + a sin(w), (110)
which is equivalent to

x = u + a cosh(v) sin(u) (111)
y = v + a sinh(v) cos(u), (112
) where the boundaries of the channel are given by v = ±V . The height of the channel is single valued for α = a cosh(V ) < 1. An example of such a channel is shown in figure 2(a).

In this case Zwanzig showed that

D e D = 1 1 + α 2 tanh(V ) V , ( 113
)
which can be easily verified using equation (97). The existence of this exact result serves as a useful verification of the general formulas derived here for the diusion equation (equations ( 24)-( 27)) along with our numerical resolution of the associated partial dierential equations. The eective diusion constant is shown in figure 2(b) for a channel with periodicity 2π with V = 0.5 for values of a taken between 0 (the flat channel) through to the maximal value of a for which the channel height is a single valued function-the agreement is seen to be perfect.

Linear channel

Here we consider dispersion in channels with piece-wise linear profiles. Such channels are of interest as they are static/non-active realizations of the saw-tooth potentials often used to study ratchets and also because of the existence of the piece-wise linear approximation scheme of [20], re derived here as equation (87). For concreteness we consider two linear channel profiles. We consider the periodic symmetric saw-tooth height profile h(x) = εζ(x) where ζ(x) = a + |x| for x ∈ [-1/2, 1/2] and shown in figure 3(a), already studied by [46][47][48][49], as well as the ratchet-like height profile ζ(x) = a + x on [0, 1] shown in figure 3(b). In the first case the channel height is continuous and equation (87) thus predicts that for the saw-tooth channel we have

D e D = 2 tan -1 (ε) ε(4a + 1) ln(1 + 1 2a ) . ( 114
)
However, for the ratchet channel the discontinuity in ζ means that equation (87) has some ambiguity. To solve the ambiguity we replace the discontinuity at x = 0, 1 by defining a finite size x 0 and taking ζ to be a straight line between (x 0 , a + x 0 ) and (1, a) (at the end of the calculation we take the limit x 0 → 1). This entails writing

ζ(x) = a + x x∈ [0, x 0 ] a + α(1 -x) x ∈ [x 0 , 1] , ( 115
)
where α = x 0 1-x 0 by continuity at x 0 . Taking the limit x 0 → 1 then recovers the ratchet channel. For a given value of x 0 we find

εζ (x) εζ(x) tan -1 (εζ (x)) s = ε tan -1 (ε) x 0 0 dx a + x + αε tan -1 (εα) 1 x 0 dx α(1 -x) + a = ε tan -1 (ε) ln 1 + x 0 a + ε tan -1 (εα) ln 1 + x 0 a , ( 116
)
and taking the limit x 0 → 1 we find,

εζ (x) εζ(x) tan -1 (εζ (x)) s = ε 2 π + 1 tan -1 (ε) ln 1 + 1 a , ( 117
)
which yields

D e D = 1 ε a + 1 2 2 π + 1 tan -1 (ε) ln 1 + 1 a . ( 118
)
The comparison of the piece-wise linear approximation with the exact results is shown in figure 3(c) we see that the approximation clearly fails for the ratchet channel, since the deviations from the basic FJ approximation are not well captured by the piece-wise linear approximation. In turn, the piece-wise linear approximation is reasonable for the saw-tooth profile for values of ε of order 1. This may linked to the fact that this approximation is exact at order ε 2 . In both cases, the piece-wise linear approximation underestimates the eective diusion constant for large ε. Indeed, in general the piece-wise approximation predicts that the eective diusion constant behaves as D e ∼ 1/ε for large ε. This is at variance with the exact numerical calculation and the large ε analysis of section 7 where the eective diusion constant saturates to a nonzero asymptotic limit (D e → a/(a + 1/4) for the saw-tooth and D e → a/(a + 1/2) for the ratchet). For the parameters a = 1/4 for the saw tooth case and a = 1/2 for the ratchet, the predicted large ε result is D e → 1/2 for both cases, in perfect agreement with the numerical results shown in figure 3.

Smooth channel

In this section we consider the case of a smoothly varying channel where all terms appearing in the perturbation series solution for the eective diusion constant D e , derived in section 4.2 are finite. In particularly we will show how a formula, which fits the numerical evaluation for the diusion constant, can be derived for all values of ε using a Padé type re-summation of the small ε perturbation expansion along with the first two terms (the O (1) and O(1/ε) ones) in the large ε expansion derived in section 7. In addition we will demonstrate numerically certain of the key results derived in this paper in the large and small ε regimes.

For the purposes of this section we consider the channel defined by

ζ(x) = 1.5 + 0.266 [cos(2πx) + sin(6πx)] . ( 119
)
In order to compare numerical results for dierent values of ε we define Y = y/ε. Shown in figures 4(a)-(c) is the numerically evaluated function f (x, Y ) for values of ε given respectively by ε = 0.01, 0.5 and 100. First, at small values of ε we see that the solution f (x, Y ) shown in figure 4(a) is independent of Y, this is the underlying hypothesis behind the FJ approximation and we see that for this narrow channel the approximation is works well. On increasing ε to 0.5, in figure 4(b), we see that the solution develops a dependence on Y. Finally for large ε, figure 4(c), we see that the analysis of section 7 is confirmed. Firstly we see that the solution vanishes within the part of the channel |Y | < ζ min to high numerical accuracy, while the solution becomes, as for the small ε case, independent of Y in the region |Y | > ζ min and takes the approximate form f (x, y) ≈ x.

We now consider the accuracy of the perturbation theory developed in this paper based on the small ε expansion, which for the smooth channel considered here is always well defined. Shown in figure 5(a) is the result of the perturbation theory starting at 24)-( 27) for the saw tooth channel for a = 0.25 (circles) and the ratchet channel for a = 0.5 (squares) compared the Kalinay-Percus piece-wise linear approximations given respectively by equations (114) (solid line) and (118) (dashed line).

order O(ε 2 ), through to that at O(ε 6 ) along with the exact numerical evaluation of the diusion constant. The results clearly improve with increasing the order of the perturbation expansion-however, they all depart strongly from the exact result when ε is O (1). Furthermore the perturbation expansion appears to give lower bounds when taken at O(ε 2 ) and O(ε 6 ) and upper bounds when taken at O(ε 4 ) (without any resummation). This alternating structure in the perturbation expansion is also seen for diusion in random Gaussian potentials [50,51]. The problem with this alternating structure of the perturbation theory is that the diusion constant will either become negative or greater than D (which is clearly not possible on physical grounds) depending on the order of perturbation theory employed. In addition, we have seen that D e saturates at large values of ε which is incompatible with a pure polynomial expansion in ε. Based on the O(ε 2 ) result of perturbation theory a number of authors have proposed the following alternative forms, which of course agree at O(ε 2 ) with the basic perturbation theory result. In our notation these alternative expressions are

D e D = 1 ζ s ζ -1 s 1 1 + ε 2 ζ 2 /ζ s 3 ζ -1 s , ( 120
)
proposed in [START_REF] Zwanzig | [END_REF]. Another alternative is given by

D e D = 1 ζ s (1 + ε 2 ζ 2 ) 1 3 ζ -1 s , (121)
which was given in [18]. A third alternative is given by equation (87), originally proposed in [20]. The comparison of these approximations with the exact numerical evaluation for the diusion constant for the case in hand is shown in figure 5(b). The advantage of these formulations is that the eective diusion constant is now bounded within the region [0, D]. However, in the limit of large ε all these proposed forms predict that the diusion constant becomes zero and all are significantly dierent from the exact result for ε of order 1. The last question we will therefore address here is whether there is a way to tie up the small and large ε results derived here to obtain a formula describing the eective diusion constant for all values of ε? Below we will pursue a Padé like approximant scheme which exploits both regimes of the perturbation theory studied in this paper. We begin by writing the small ε expansion for the eective diusivity given by equation ( 67)

D e D = 1 ζ s ζ -1 s 1 - ζ 2 /ζ s 3 ζ -1 s ε 2 + ζ 2 /ζ 2 s 9 ζ -1 2 s + 4 ζ 4 /ζ s + ζ 2 ζ s 45 ζ -1 s ε 4 + ... = D FJ 1 + λ 2 ε 2 + λ 4 ε 4 + ... , ( 122
)
where λ i = C i C 0 . In principle one can gain in precision by also using the O(ε 6 ) correction, though the resulting algebraic complication does not really justify the gain in accuracy. Using the large ε expansion developed in section 7 we also write A Padé type re-summation is constructed by writing

D e D = D FJ 1 + a 1 ε + a 2 ε 2 + a 3 ε 3 1 + b 1 ε + b 2 ε 2 + b 3 ε 3 , ( 124
)
where the a i and b i are obtained from the two expansions, for small and large ε. The coecients are found to be given by (

) 125
The results of this approximation for the class of channels described by equation (119) is shown in figure 5(b). We see that the resulting approximation is in very good agreement with the exact numerical results and that, as is to be expected, the maximum errors arise in the crossover region where ε is of order 1.

Conclusion

In this paper we have examined the eective dispersion of Brownian particles in a channel geometry. Rather than attempting to write an eective one-dimensional diusion equation to describe the system as in many past studies, we have exploited a direct formula for the eective diusion constant which is obtained from the system of time independent equations equations (24)-( 27) which were derived here for completeness but are in agreement with a general theory of dispersion coecients developed in [34] and [35] and also with the macro transport theory of Brenner [START_REF] Brenner | Macrotransport Theory[END_REF] when applied to the channel problem [28]. Using a formulation in terms of analytic functions (equations ( 24)-( 27)), we were able to develop a compact formulation of the perturbation theory in channel width scale ε. Notably we were able to show that the homogenization scheme proposed by Kalinay and Percus [20] gives the correct diusion constant at order ε 6 . The exact results of Zwanzig for parameterized channels and the piece-wise linear approximation of [20] can be deduced from the formalism developed here in a very compact manner. Most studies on diusion in channels are focused on narrow channels. The formalism developed here was also exploited to study wide channels, yielding a physically intuitive result in the limit or very large channels, but also giving the first non-trivial correction which takes a remarkably universal form. Finally we were able to show how the small and large ε expansions for the eective diusion constant can be combined in a Padé type approximant which accurately describe continuous channels for all values of ε.

There are a number of potential extensions of this work based on the static approach proposed in [34,35] which can be extended to systems having a current in their equilibrium states. This means that the approach could be applied to dispersion in channels with driving forces both along the channel and normal to the channel [11,52,53] (arising from gravity for example). Furthermore the approach of [34,35] can also be adapted to extract the finite time behavior of transport coecients and could be applied to the problem of diusion in channels.

Figure 1 .

 1 Figure 1. (a) The domain of the analytic function w(z). (b) The domain of the analytic function w(z) = w(z) -z. The wiggly lines represents a jump of +1 in the solution. (c) The domain of the analytic function w(ξ), where ξ = exp [-iπ(z -ih min )].A jump of +1 in the solution occurs between the two curved lines.

Figure 2 .

 2 Figure 2. (a) Example of a channel in the class parameterized by equations (111) and (112), where a = 0.5 and V = 0.5. (b) Comparison of the numerical evaluation of the system of equations (24)-(27) (points) for the eective diusion constant of the channel parameterized by equations (111) and (112) and the exact result equation (113) (solid line).

Figure 3 .

 3 Figure 3. (a) Schematic of continuous symmetric piece-wise linear channel (sawtooth shape) (b) Schematic of discontinuous ratchet-like channel. (c) Comparison of the (exact) numerical evaluation of the system of equations (24)-(27) for the saw tooth channel for a = 0.25 (circles) and the ratchet channel for a = 0.5 (squares) compared the Kalinay-Percus piece-wise linear approximations given respectively by equations (114) (solid line) and (118) (dashed line).

  min + ln 2 πε + ... = D ∞ 1 + α ε + ... . (123)

Figure 4 .

 4 Figure 4. (a) Numerical solution of f (x, Y ) for ε = 0.01. The function is independent of the coordinate Y = y/ε in accordance with the FJ approximation. The resulting diusion constant is thus well approximated by D e ∼ D FJ given in equation (47). (b) Numerical solution of f (x, Y ) for ε = 0.5. The function is no longer independent of Y and is also non-zero in the region of the channel |Y | < ζ min . (c) Numerical solution of f (x, Y ) for ε = 100. In agreement with the results of section 7 we see that f = 0 in the region of the channel |Y | ζ min and is independent of Y for |Y | > ζ min .

Figure 5 .a 1 = b 1 = [λ 2 2 + 2 a 2 = λ 2 - λ 4 λ 2 b 2 = - λ 4 λ 2 a 3 = 2 + λ 4 (K - 1 )

 5112222223241 Figure 5. (a) Exact numerical value of D e (circles) for the class of channels described by equation (119) as a function of ε compared with perturbation theory in ε 2 at order 1 (dot-dashed line), 2 (dashed line) and 3 (solid line). (b) Same exact numerical values for D e as in (a) (circles) compared with the O(ε 2) re-summations given by equation (120) due to Zwanzig[START_REF] Zwanzig | [END_REF] equation (121) (dashed line), due to Reguera and Rubi[18] (dot-dashed line) and equation (87) due to Kalinay and Percus[20] (dotted line). Shown as a solid line is the Padé approximant scheme proposed here in the system of equations (124)-(125).

Appendix. Derivation of equation (41)

The function w(x + iy) is analytic, thus for z = x + iy, we have for all closed paths C C dzw(z) = 0.

(A.1)

Note that here we can add an arbitrary constant to w(z) and of course the above remains trivially true (justifying making the choice v 0 = 0 in equation ( 36)). If we choose C to be the boundary of one periodic cell, with x ∈ [0, L], we get

idy[u(L, y) + iv(L, y)]

idy[u(0, y) + iv(0, y)].

(A.2) If we use the symmetries u(x, -y) = u(x, y) and v(x, -y) = -v(x, y), and the periodicity in x of all functions, the expression becomes

From equation (40), we obtain

and the boundary condition (38) yields .5)