THE JOURNAL OF CHEMICAL PHYSICS 142, 014106 (2015)

Generalized method calculating the effective diffusion coefficient in periodic channels Pavol Kalinaya) Institute of Physics, Slovak Academy of Sciences, Dúbravská cesta 9, 84511 Bratislava, Slovakia

(Received 29 October 2014; accepted 15 December 2014; published online 7 January 2015) The method calculating the effective diffusion coefficient in an arbitrary periodic two-dimensional channel, presented in our previous paper [P. Kalinay, J. Chem. Phys. 141, 144101 (2014)], is generalized to 3D channels of cylindrical symmetry, as well as to 2D or 3D channels with particles driven by a constant longitudinal external driving force. The next possible extensions are also indicated. The former calculation was based on calculus in the complex plane, suitable for the stationary diffusion in 2D domains. The method is reformulated here using standard tools of functional analysis, enabling the generalization. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4905079]

I. INTRODUCTION

net flux

Study of transport through nonhomogeneous quasi onedimensional (1D) systems, like pores, channels in nanomaterials or biological structures, is often considerably facilitated by use of an effective 1D description. Typically, the motion of particles along the channel is of our primary interest. Then, it is natural to look for equations describing evolution of the system only in the longitudinal coordinate x. A prototype of such equations is the generalized FickJacobs (FJ) equation1–3 ∂t p(x, t) = ∂x A(x)D(x)∂x [p(x, t)/A(x)],

(1.1)

governing evolution of the 1D (or linear) density of particles p(x, t). The functions A(x) and D(x) involve properties of the channel in the transverse directions y, which cannot be neglected, as they significantly influence the longitudinal transport. In general, A(x) is the local configuration integral (at a fixed x) and D(x) is an effective diffusion coefficient. In the simplest case, diffusion in a 2D channel bounded by a function y = h(x) > 0 and the x axis, A(x) = h(x) is width of the channel. In equilibrium, the 2D density ρ(x, y, t) is constant, hence the 1D density, defined as  h(x) p(x, t) = ρ(x, y, t)dy, (1.2) 0

is proportional to A(x) and so Eq. (1.1) gives correctly ∂t p(x, t) = 0. Introducing A(x) in Eq. (1.1) maintains the 1D mass conservation along the x axis in general. If A(x) is interpreted as Boltzmann’s factor of an effective “entropic” potential Ue f f (x), A(x) = exp[− βUe f f (x)], β = 1/k BT denotes the inverse temperature, Eq. (1.1) for a constant D(x) becomes the Smoluchowski equation. The effective diffusion coefficient D(x) reflects variations of the longitudinal diffusion rate due to confinement.2,3 It is thoroughly determined by the stationary flow,4,5 where the a)Electronic mail: [email protected]

0021-9606/2015/142(1)/014106/6/$30.00

 J(x, t) =

h(x)

j x (x, y, t)d y

(1.3)

0

⃗ is constant, J(x, t) = J , 0; ⃗j(x, y, t) = −D0∇ρ(x, y, t) denotes the 2D flux density and D0 is the bulk diffusion constant. Equation (1.1) represents the 1D mass conservation, ∂t p(x, t) = −∂x J(x, t), hence J(x, t) = −A(x)D(x)∂x [p(x, t)/A(x)].

(1.4)

The stationary density ρ(x, y) is given by solution of the 2D diffusion equation ( ) ∂t ρ(x, y) = D0 ∂x2 + ∂y2 ρ(x, y) = 0, (1.5) supplemented by the reflecting (Neumann) boundary conditions (BC) ∂y ρ(x, y)| y=h(x) = h ′(x)∂x ρ(x, y)| y=h(x), ∂y ρ(x, y)| y=0 = 0,

(1.6)

up to two integration constants, J and (an irrelevant) equilibrium density ρ0. Thus, if the solution ρ(x, y) is known, the corresponding p(x) and J can be calculated according to Eqs. (1.2) and (1.3), respectively, and D(x) is finally expressed from Eq. (1.4). The algorithms calculating D(x) are based mainly on scaling of the transverse lengths y, h(x) by a small parameter √ ϵ → 0, which is equivalent to introducing anisotropy of the diffusion constant in Eq. (1.5); its transverse value is D0/ϵ → ∞. If the transverse relaxation was infinitely fast, ρ(x, y) would become constant (equilibrated) across the channel, so ρ(x, y) = p(x)/A(x) to satisfy Eq. (1.2). Then j x (x, y) = −D0∂x [p(x)/A(x)], applied in Eqs. (1.3) and (1.4), would result in D(x) = D0. In real channels, the transverse diffusion is not so fast and ρ(x, y) deviates from p(x)/A(x), causing D(x) < D0. Series of corrections to D0 controlled by ϵ can be obtained by recurrence schemes;6–9 the result for the considered

142, 014106-1

© 2015 AIP Publishing LLC

014106-2

Pavol Kalinay

J. Chem. Phys. 142, 014106 (2015)

2D channel is  ϵ ϵ D(x)/D0 = 1 − h ′2(x) + h ′(x) 9h ′3(x) 3 45  +h(x)h ′(x)h ′′(x) − h2(x)h(3)(x) + ··· √ √ (1.7) ≃ arctan[ ϵ h ′(x)]/[ ϵ h ′(x)] 2

after neglecting h ′′(x) and the higher derivatives and summing the series to infinity. Validity of the generalized FJ equation (1.1) with D(x) calculated according to the formula (1.7) was verified up to |h ′(x)| ∼ 1 by the Brownian simulations.10–12 Still, its application becomes questionable especially in geometries with cusps13,14 or jumps,15–19 where the derivatives contained in Eq. (1.7) cannot be expressed. A recent work20 offers an alternative way for calculation of D(x) in periodic 2D channels. It is based on solution of the Laplace equation (1.5) in the complex plane. The stationary density ρ(x, y) is identified with the imaginary part of an analytic complex function f (z) = ibz + 2i

∞ 

[an sin(nk z) + bn cos(nk z)],

(1.8)

n=1

k = 2π/L, z = x + i y, and L is period of the channel. The coefficients an , bn , and b are calculated from the condition ∞    1 = −bh(x) + 2 bn sin(nk x) − an cos(nk x)

of the stream lines, Fig. 1, i.e., the flux density ⃗j(x, y) = ⃗j(x + L, y). For diffusion, we have jα (x, y) = −D0∂α ρ(x, y) = −D0∂α ρ(x + L, y),

α = x, y,.... For the longitudinal component, α = x, the condition can be satisfied only for ∞   ρ(x, y) = ρ( ¯ y)x + ρ¯ 0( y) + ρ n ( y)sin(nk x) n=1

 + ρ¯ n ( y)cos(nk x) ,

n=1

 + ρ¯ n ( y)cos(nk x) .

(2.3)

The solution of the zeroth term ρ¯ 0( y) = a + a¯ y, the particular solutions of [∂y2 −(nk)2]ψ( y) = 0 are sinh(nk y) and cosh(nk y). BC (1.6) at y = 0 excludes sinh(nk y) and set a¯ = 0, so the 2D density has the form ∞  ρ(x, y) = a+bx +2 (an sinnk x +bn cosnk x)cosh(nk y); n=1

(2.4) (1.9)

by the Fourier transform (FT). They are finally applied in the formula for D(x), derived using Eq. (1.4), ∞   sinh[nk h(x)] D0 = −h(x)∂x bx + 2 D(x) nk h(x) n=1  × (an sinnk x + bn cosnk x) ,

(2.2)

k = 2π/L and the functions ρ n ( y), ρ¯ n ( y) ρ( ¯ y) are fixed from the next conditions. The transverse component of Eq. (2.1) requires L∂y ρ( ¯ y) = 0, hence ρ( ¯ y) = b constant. Applying the diffusion equation (1.5) on the stationary ρ(x, y), Eq. (2.2), we obtain ∞    0 = ∂y2 ρ¯ 0( y) + ∂y2 − (nk)2 ρ n ( y)sin(nk x)

n=1

× sinh[nk h(x)]

(2.1)

(1.10)

corresponding to Eqs. (2.13)–(2.15) in Ref. 20. The derivatives of h(x) do not enter the procedure in any step, so the method is also suitable for nonanalytic periodic functions h(x), i.e., for the channels with cusps or jumps. The problem appears in the formulation, which uses properties of the complex analytic functions, so it is applicable only for 2D diffusion. The present paper reformulates the method in a way enabling its extension for 3D channels, which is done in Sec. II. Section III generalizes it for diffusion of particles driven by a constant longitudinal force F. Another possible extensions of the method are discussed in the Conclusion.

we retain notation from Ref. 20. Notice that this formula corresponds to Im[ f (z)], Eq. (1.8), up to the constant a, representing an equilibrium “offset” of the density, which does not contribute to the net flux. Indeed, applying (2.4) in Eqs. (1.3) and (1.2), we find ∞    J/D0 = −bh(x) + 2 bn sin(nk x) − an cos(nk x) n=1

× sinh[nk h(x)],

(2.5)

∞  p(x) = (a + bx)h(x) + 2 (an sinnk x + bn cosnk x) n=1

× sinh[nk h(x)]/(nk).

(2.6)

According to Eq. (2.5), the coefficients an , bn , and b are proportional to J/D0. Setting J/D0 = 1 and applying Eq. (2.6) in

II. REFORMULATION OF THE METHOD

We reformulate here the procedure calculating D(x), indicated in the Introduction, without use of the complex analytic function f (z), Eq. (1.8) and verify the formulas (1.9) and (1.10). Next, the method is extended to 3D geometries with cylindrical symmetry. The key symmetry determining the stationary flow in a periodic channel of period L, h(x) = h(x + L), is periodicity

FIG. 1. A periodic channel defined by y = h(x) = h(x + L). Thin lines describe the stream lines, which are also periodic, so the flux density ⃗j(x, y) = ⃗j(x + L, y).

014106-3

Pavol Kalinay

J. Chem. Phys. 142, 014106 (2015)

calculation of D(x) from Eq. (1.4), we restore the formula (1.10) and Eq. (2.5) becomes identical to the condition (1.9). The condition (1.9) is equivalent to BC (1.6) at y = h(x). Substituting Eq. (2.4) for ρ(x, y) in (1.6), we get 0=

∞ 

2nk(an sinnk x + bn cosnk x)sinh[nk h(x)] − h ′(x)

n=1

   ∞ × b+ 2nk(an cosnk x − bn sinnk x)cosh[nk h(x)] .  n=1  (2.7) The same equation is obtained from the derivative of Eq. (1.9) according to x. Of course, the condition (1.9) does not contain h ′(x) and so it is preferable for calculation of D(x) in the channels with cusps or jumps. Extension to 3D channels with cylindrical symmetry is now straightforward. The stationary 3D density ρ(x,r) does not depend on the angle φ in the cylindrical coordinates (x,r,φ) and obeys the diffusion equation ) ( 2 1 (2.8) ∂t ρ(x,r) = D0 ∂x + ∂r r∂r ρ(x,r) = 0 r with the reflecting BC ∂r ρ(x,r)|r =h(x) = h ′(x)∂x ρ(x,r)|r =h(x);

(2.9)

h(x) denotes now the local radius of the channel. Periodicity of the x component of ⃗j, j x (x,r) = −D0∂x ρ(x,r) determines ρ(x,r) of the form (2.2) with y replaced by r; periodicity of the radial component jr (x,r) = −D0∂r ρ(x,r) fixes ρ(r) ¯ = b. Application of Laplacian on ρ(x,r) according to Eq. (2.8) gives ∞    1 1 ∂r r∂r ρ¯ 0(r) + ∂r r∂r − (nk)2 ρ n (r)sin(nk x) r r n=1  + ρ¯ n (r)cos(nk x) . (2.10)

applied on Eq. (2.11) results in p(x, t) = (a + bx)A(x) + 2π

∞  h(x) n=1

nk

I1[nk h(x)]

× (an sinnk x + bn cosnk x).

(2.15)

According to Eq. (2.13), the coefficients an , bn , and b are proportional to J/D0. If we set J/D0 = 1, the final formula for D(x), calculated from Eq. (1.4), reads ∞   D0 I1[nk h(x)] = −πh2(x)∂x bx + 2 D(x) nk h(x) n=1  × (an sinnk x + bn cosnk x) .

(2.16)

For a flat channel, h(x) = h constant, Eq. (2.13) gives b = −1/(πh2), an ,bn = 0. So D(x)/D0 = 1 from Eq. (2.16), as expected. One can easily check again that the derivative of Eq. (2.13) according to x is equivalent to the BC (2.9) for ρ(x,r) defined by Eq. (2.11). The Eqs. (1.10) and (2.16) for D(x) simplify the LifsonJackson formula21 for the effective diffusion coefficient De f f averaged over one period of the channel,  L  L 1 D0dx D0 = 2 A(x)dx De f f L 0 0 A(x)D(x)  L  L 1 p(x) =− 2 A(x)dx A(x) 0 L 0  L b =− A(x)dx. (2.17) L 0 The oscillating terms in the sums subtract one another at the limits x = 0 and L and so De f f is determined solely by the coefficient b.

0=

The zeroth term is solved by ρ(r) ¯ = a + alnr; ¯ the particular solutions of ρ n (r), ρ¯ n (r) are the modified Bessel functions I0(nkr) and K0(nkr). The 3D density has no singularity at the x axis, hence a¯ = 0 and also K0(nkr) is excluded. Then, ρ(x,r) is of the form ρ(x,r) = a + bx +

∞  (an sinnk x + bn cosnk x)I0(nkr).

(2.11)

n=1

The net flux in 3D symmetric channels,  h(x) J(x, t) = 2πr j x (x,r,t)dr,

(2.12)

0

calculated for ρ(x,r), Eq. (2.11), gives the condition J/D0 = −bA(x) + 2π

∞  

bn sin(nk x) − an cos(nk x)

A(x) = πh2(x) here. The 1D density,  h(x) p(x, t) = 2πr ρ(x,r,t)dr, 0

The algorithm presented in Sec. II is also applicable for diffusion in a time independent field, which does not violate periodicity of the stationary flux density, ⃗j(x, y) = ⃗j(x + L, y). In this section, we demonstrate calculation of D(x) for diffusion in 2D and 3D channels driven by a constant longitudinal force F. The stationary 2D density ρ(x, y) obeys now the Smoluchowski equation,   ∂t ρ(x, y) = D0 ∂x e−βU (x)∂x e βU(x) + ∂y2 ρ(x, y) = 0; (3.1) β = 1/k BT denotes the inverse temperature and the potential U(x) = −F x generates the constant longitudinal force F. The flux density is then given by
⃗j(x, y) = −D0 e β F x ∂x e−β F x ρ(x, y),∂y ρ(x, y) , (3.2)



n=1

× h(x)I1[nk h(x)];

III. DIFFUSION DRIVEN BY A LONGITUDINAL FORCE

(2.13)

as Eq. (3.1) represents the mass conservation. It is supplemented by the reflecting BC ∂y ρ(x, y)| y=h(x) = h ′(x)e β F x ∂x e−β F x ρ(x, y)| y=h(x), ∂y ρ(x, y)| y=0 = 0

(2.14)

at the upper and lower boundaries.

(3.3)

014106-4

Pavol Kalinay

J. Chem. Phys. 142, 014106 (2015)

Periodicity of the x-component j x (x, y) requires ρ(x, y) = ρ( ¯ y)e β F x +

∞ 

ρ n ( y)eink x ;

(3.4)

n=−∞

then we get j x (x, y) = D0 n ( βF − ink)ρ n ( y)eink x = j x (x + L,y) from Eq. (3.2), k = 2π/L. Periodicity of j y (x, y) fixes ρ( ¯ y) = a constant. The equations for ρ n ( y) are obtained after substituting ρ(x,y) (3.4) in Eq. (3.1), 

0=

∞  

 ∂y2 − (nk)2 −ink βF ρ n ( y)eink x .

(3.5)

n=−∞

¯ and the particular soluThen the zeroth term ρ0( y) = b + by tions for n , 0 are exp[± (nk)2 +ink βF y]. The BC (3.3) at y = 0 sets b¯ = 0 and allows only cosh combination of both exponentials, resulting in   ρ(x, y) = ae β F x + b+ cn cosh[ (nk)2 +ink βF y]eink x .

density ρ(x,r) obeys ) ( βF x −β F x 1 ∂t ρ(x,r) = D0 ∂x e ∂x e + ∂r r∂r ρ(x,r) = 0 r

(3.9)

and Neumann’s BC. Periodicity of ⃗j(x,r) requires ρ(x,r) of the form ∞  ρ(x,r) = ae β F x + ρ n (r)eink x , (3.10) n=−∞

the functions ρ n (r) are found from the Smoluchowski Eq. (3.9),  ∞   1 2 ∂r r∂r − (nk) −ink βF ρ n (r)eik r . (3.11) 0= r n=−∞ ¯ The zeroth component is solved by ρ0(r) = b+ blnr, the particular solutions for  n , 0 are the modified Bessel functions I0(z), K0(z) of z = (nk)2 +ink βFr. The 3D density has no singularity at r = 0, hence    ρ(x,r) = ae β F x + b+ cn I0 (nk)2 +ink βFr eink x , n,0

n,0

(3.6)

(3.12)

The coefficients  cn are now complex. We can write cn = bn −ian and denote (nk)2 +ink βF = κ n +iη n , κ −n = κ n , and η −n = −η n . Requirement of real ρ(x, y) confines b−n = bn and a−n = −an . Using an , bn , the sum in Eq. (3.6) approaches the sum in Eq. (2.4) for F → 0. The net flux is expressed using Eqs. (1.3) and (3.2),

cn are complex constants. Then, the net flux J is calculated according to Eq. (2.12) with j x (x,r) the same as in Eq. (3.2),  ( βF −ink) J/D0 = πh2(x) βFb+ 2π cn  2 +ink βF (nk) n,0  × h(x)I1[ (nk)2 +ink βF h(x)]eink x . (3.13)

J/D0 = βFbh(x) +



cn 

( βF −ink)

(nk)2 +ink βF n,0   × sinh (nk)2 +ink βF h(x) eink x .

(3.7)

Its derivative according to x is identical to the BC (3.2) for ρ(x,y) (3.6) at y = h(x). The coefficients b and cn are fixed unambiguously from FT of this equation, while a represents an equilibrium “offset,” which does not contribute neither to the flux, nor to D(x). After integration of p(x) according to Eq. (1.2), the final formula for D(x) is obtained from Eq. (1.4). Here, the configuration integral A(x) contains the Boltzmann factor, A(x) = h(x)e β F x , hence   D0 = −h(x)e β F x ∂x e−β F x p(x)/h(x) D(x)   = −h(x)e β F x ∂x be−β F x + cn eink x−β F x n,0

  sinh (nk)2 +ink βF h(x)  ×  ; (nk)2 +ink βF h(x)

(3.8)

we set J/D0 = 1 here as well as in Eq. (3.7), fixing b and cn . In a flat channel, h(x) = h constant, βFb = 1/h, and cn = 0; p(x) = [ae β F x + 1/( βF h)]. Eq. (3.8) then gives D(x) = D0, as expected. For completeness, we add briefly the calculation of D(x) for driven diffusion in 3D channels of cylindrical symmetry, bounded by radius r = h(x). The symmetric 3D stationary

Finally, p(x) is calculated using Eq. (2.14) and D(x) is expressed from Eq. (1.4). Again, A(x) contains the Boltzmann weight, A(x) = πh2(x)e β F x , so   D0 = −πh2(x)e β F x ∂x be−β F x + 2 cn eink x−β F x D(x) n,0   2 I1[ (nk) +ink βF h(x)] ×  ; (3.14) (nk)2 +ink βF h(x) the coefficients b, cn are fixed for J/D0 = 1. Derivation of the effective mobility µ∗ and diffusion coefficient D ∗ averaged over one period requires us to apply the treatment of Reimann et al.,22,23 modified for the generalized FJ equation (1.1) with D(x) , 1. The mean velocity ⟨ x⟩ ˙ = µ∗F ∗ and D are expressed by the formulas ⟨ x⟩ ˙ =

L , T1(x 0 → x 0 + L)

D∗ =

L 2 ∆T2(x 0 → x 0 + L) , 2 [T1(x 0 → x 0 + L)]3

(3.15)

using the moments Tl (x 0 → ξ) of the time in which a particle inserted at x 0 reaches the position ξ,  ∞  ξ Tl (x 0 → ξ) = − t l dt∂t p(x,t|x 0,0)dx 

0 ∞

=l



−∞ ξ

p(x,t|x 0,0)dx.

t l−1dt 0

(3.16)

−∞

We suppose F > 0, x 0 < ξ; the conditional probability p(x,t| x 0,0) of finding a particle at (x, t) if it started from x 0 obeys

014106-5

Pavol Kalinay

J. Chem. Phys. 142, 014106 (2015)

Eq. (1.1). The corresponding solution of the problem, satisfying Tl (ξ → ξ) = 0, is  ξ  x dx Tl (x 0 → ξ) = l A( y)Tl−1( y → ξ)d y, x 0 A(x)D(x) −∞ (3.17) ∆T2 = T2 −T12

T0(x 0 → ξ) = 1 for any x 0 < ξ. We use l = 1,2 and in Eqs. (3.15). Let us introduce the cross section area σ(x) = h(x) or πh2(x) in 2D or 3D channels, respectively, then A(x) = σ(x) e β F x . For any dimensionality, D(x) is of the form   D0/D(x) = −σ(x)e β F x ∂x e−β F x Q(x) ,

(3.18)

see Eqs. (3.8) and (3.14). The function Q(x) = Q(x + L) is periodic with the constant (leading) term b. Then Tl and so ⟨ x⟩ ˙ and D ∗ can be expressed in the following way:  ξ   1 dx∂x e−β F x Q(x) T1(x 0 → ξ) = − D0 x 0  x × e β F y σ( y)d y −∞



1 = D0

L

σ(x)Q(x)dy

(3.19)

0

for ξ = x 0 + L according to Eq. (3.17), integrating by parts and using periodicity of σ(x), Q(x). T1(x 0 → x 0 + L) = T1(0 → L) does not depend on x 0. x first formula, −∞ e β F y σ( y)dy = (1 − e−β F L )−1  xIn the × x−L e β F y σ( y)dy.22 For F → 0, the integral over one period L becomes 0 σ( y)d y, independent of x. The double integral decouples and if applied in Eq. (3.15), it results in ⟨ x⟩ ˙ = βF D LJ with Lifson-Jackson’s diffusion coefficient D LJ restored, taking 1 − e−β F L ≃ βF L. In general, one obtains ( ⟨ x⟩ ˙ = D0 L

) −1

L

σ(x)Q(x)dx

.

(3.20)

0

Its proportionality to F for small F is determined by Q(x); e.g., in a flat channel with σ(x) = σ constant, Eq. (3.7) or (3.13) for J = D0 give b = 1/( βFσ). Applied in Eq. (3.20), we find ⟨ x⟩ ˙ = βF D0. A similar algebra is used to obtain D ∗ from Eq. (3.15). Following Reimann’s Appendix A22 and applying Eq. (3.18), we get  ξ   2 ∆T2(x 0 → ξ) = dx∂x e−β F x Q(x) 3 D02(1 − e−β F L ) x0  z 2  x   × dz∂z e−β F z Q(z) e β F y σ( y)dy x−L z−L  L 2 = 2 dxσ(x)Q2(x) D0 (1 − e−β F L ) 0  x −β F x ×e e β F y σ( y)dy (3.21) x−L

for ξ = x 0 + L. The quadruple integral decouples in the limit F → 0, giving Lifson-Jackson’s formula if applied in Eq. (3.15).

For a general F, we obtain D = ∗

L 2 D0

L 0

x dxσ(x)Q2(x)e−β F x x−L e β F y σ( y)d y . 3  L (1 − e−β F L ) 0 dxσ(x)Q(x)

(3.22)

Knowing Q(x), i.e., the coefficients b and cn , solving Eq. (3.7) or (3.13) with J/D0 = 1, one should find the effective µ∗ and D ∗ from Eqs. (3.20) and (3.22) exactly in principle for any shape of the periodic channel and the driving force F. Of course, the algorithm presented here is suitable mainly for numerical solution; the results will be influenced by the size (truncation) of the used basis and numerical errors. Still, construction of a compact formula approximating Q(x) based on the presented analysis, as well as its comparison to the expansions of µ∗, D ∗ in the scaling parameter ϵ 24 remain a challenge for our next study.

IV. CONCLUSION

The aim of this paper was to generalize the method of calculation of the effective diffusion coefficient D(x) in periodic channels, suggested in Ref. 20. In contrast to the widely used formulas derived by scaling of the transverse lengths and homogenization schemes,3,6–9 it does not need to express and use derivatives of the function h(x) shaping the channel. So it is also applicable for the channels with cusps or jumps of their width. The method is based on a close connection of D(x) to the stationary flow. If the stationary density ρ(x, y, t) = ρ(x, y) is known, D(x) is fixed unambiguously applying the relations (1.2)–(1.4). The previous formulation20 expressed ρ using an analytic complex function f (x + i y), Eq. (1.8), where ρ(x, y) = Im[ f (x +i y)]. Of course, it is applicable only for 2D diffusion described in the complex plane. We show here that the calculation can be reformulated by standard tools of the functional analysis and so generalized for diffusion in 3D channels, as well as diffusion driven by a constant longitudinal force F. The key point of the calculation is periodicity of the stationary flux density ⃗j(x, y) = ⃗j(x + L, y), reproducing periodicity of the channel. Then, ρ(x, y) and so D(x) can be expressed in Fourier series, Eq. (1.10); the coefficients b, an , and bn are fixed unambiguously by FT of the equation for the net flux, Eq. (1.9). We derived here similar equations for 3D channels with cylindrical symmetry, as well as for the driven diffusion. We suppose that the method can be also successfully applied on systems with other potentials, which do not violate periodicity of the flux density ⃗j(x, y). An example is diffusion biased by a transverse constant (gravitational) force,25–28 or in a channel with soft walls, shaped by a periodic potential U(x, y) = U(x + L, y), e.g., a parabolic well U(x, y) = α(x) y 2 with periodic stiffness α(x) = α(x + L). The homogenization mapping procedure applied on such systems29,30 generates the corresponding effective 1D equations reducible to the form (1.1) in the limit of stationary flow. The method as presented here is appropriate mainly for numerical calculation of D(x), complementary to Brownian simulations used for studying periodic systems till now. Nevertheless, this scheme might serve as a basis for constructing

014106-6

Pavol Kalinay

an analytic integral formula for D(x), complementary to the formulas3,7,8 using h ′(x), which could be applicable also for channels with steep changes of width, h ′(x) > 1. Finding such a formula will be the object of our study in the future.

ACKNOWLEDGMENTS

Support from VEGA Grant No. 2/0015/15 is gratefully acknowledged. 1M.

H. Jacobs, Diffusion Processess (Springer, New York, 1967). Zwanzig, J. Phys. Chem. 96, 3926 (1992). 3D. Reguera and J. M. Rubí, Phys. Rev. E 64, 061106 (2001). 4P. Kalinay and J. K. Percus, Phys. Rev. E 78, 021103 (2008). 5P. Kalinay, EPJ-ST 223, 3027 (2014). 6P. Kalinay and J. K. Percus, J. Chem. Phys. 122, 204701 (2005); J. Stat. Phys. 123, 1059 (2006). 7P. Kalinay and J. K. Percus, Phys. Rev. E 74, 041203 (2006). 8S. Martens, G. Schmid, L. Schimansky-Geier, and P. Hänggi, Phys. Rev. E 83, 051135 (2011); Chaos 21, 047518 (2011). 9K. D. Dorfman and E. Yariv, J. Chem. Phys. 141, 044118 (2014). 10A. M. Berezhkovskii, M. A. Pustovoit, and S. M. Bezrukov, J. Chem. Phys. 126, 134706 (2007). 11A. M. Berezhkovskii, M. A. Pustovoit, and S. M. Bezrukov, Phys. Rev. E 80, 020904(R) (2009); Chem. Phys. 375, 523 (2010). 12I. Pineda, M. V. Vazquez, A. M. Berezhkovskii, and L. Dagdug, J. Chem. Phys. 135, 224101 (2011); 136, 204106 (2012). 13A. M. Berezhkovskii, L. Dagdug, Y. A. Makhnovskii, and V. Yu. Zitserman, J. Chem. Phys. 132, 221104 (2010); M. V. Vázquez, A. M. Berezhkovskii, 2R.

J. Chem. Phys. 142, 014106 (2015) and L. Dagdug, ibid. 129, 046101 (2008); M. V. Vázquez and L. Dagdug, J. Mod. Phys. 2, 284 (2011). 14L. Dagdug, M.-V. Vazquez, A. M. Berezhkovskii, and S. M. Bezrukov, J. Chem. Phys. 133, 034707 (2010). 15A. M. Berezhkovskii, A. V. Barzykin, and V. Yu. Zitserman, J. Chem. Phys. 131, 224110 (2009). 16F. Marchesoni, J. Chem. Phys. 132, 166101 (2010). 17L. Dagdug, A. M. Berezhkovskii, Y. A. Makhnovskii et al., J. Chem. Phys. 134, 101102 (2011); 136, 214110 (2012). 18A. E. Antipov, A. V. Barzykin, A. M. Berezhkovskii et al., Phys. Rev. E 88, 054101 (2013). 19P. Kalinay and J. K. Percus, Phys. Rev. E 82, 031143 (2010). 20P. Kalinay, J. Chem. Phys. 141, 144101 (2014); Erratum 141, 169902 (2014) 21S. Lifson and J. L. Jackson, J. Chem. Phys. 36, 2410 (1962). 22P. Reimann, C. Van den Broeck, H. Linke, P. Hänggi, J. M. Rubì, and A. Pérez-Madrid, Phys. Rev. E 65, 031104 (2002); Phys. Rev. Lett. 87, 010602 (2001). 23D. Reguera, G. Schmid, P. S. Burada, J. M. Rubì, P. Reimann, and P. Hänggi, Phys. Rev. Lett. 96, 130609 (2006). 24N. Laachi, M. Kenward, E. Yariv, and K. D. Dorfman, EPL 80, 50009 (2007). 25F. Marchesoni and S. Savel’ev, Phys. Rev. E 80, 011120 (2009); 73, 021102 (2006). 26P. S. Burada and G. Schmid, Phys. Rev. E 82, 051128 (2010). 27P. S. Burada, G. Schmid, D. Reguera, M. H. Vainstein, J. M. Rubì, and P. Hänggi, Phys. Rev. Lett. 101, 130602 (2008); Eur. Phys. J. B 69, 11 (2009). 28D. Mondal, M. Das, and D. S. Ray, J. Chem. Phys. 132, 224102 (2010); D. Mondal and D. S. Ray, Phys. Rev. E 82, 032103 (2010). 29P. Kalinay, Phys. Rev. E 80, 031106 (2009); 84, 011118 (2011). 30P. Kalinay and J. K. Percus, Phys. Rev. E 83, 031109 (2011).

Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.