Bulk-mediated diffusion on a planar surface: full solution

# Bulk-mediated diffusion on a planar surface: full solution

Aleksei V. Chechkin Institute for Theoretical Physics NSC KIPT, Akademicheskaya st.1, 61108 Kharkov, Ukraine Institute for Physics & Astronomy, University of Potsdam, D-14476 Potsdam-Golm, Germany    Irwin M. Zaid Rudolf Peierls Centre for Theoretical Physics, University of Oxford, 1 Keble Road, Oxford OX1 3NP, United Kingdom    Michael A. Lomholt MEMPHYS - Center for Biomembrane Physics, Department of Physics, Chemistry, and Pharmacy, University of Southern Denmark, Campusvej 55, 5230 Odense M, Denmark    Igor M. Sokolov Institut für Physik, Humboldt Universität zu Berlin, Newtonstraße 15, 12489 Berlin, FRG    Ralf Metzler Institute for Physics & Astronomy, University of Potsdam, D-14476 Potsdam-Golm, Germany Department of Physics, Technical University of Tampere, FI-33101 Tampere, Finland
5th July 2019
###### Abstract

We consider the effective surface motion of a particle that intermittently unbinds from a planar surface and performs bulk excursions. Based on a random walk approach we derive the diffusion equations for surface and bulk diffusion including the surface-bulk coupling. From these exact dynamic equations we analytically obtain the propagator of the effective surface motion. This approach allows us to deduce a superdiffusive, Cauchy-type behavior on the surface, together with exact cutoffs limiting the Cauchy form. Moreover we study the long-time dynamics for the surface motion.

###### pacs:
05.40.Fb,02.50.Ey,82.20.-w,87.16.-b

## I Introduction

Interfaces and the interaction of particles with them play a crucial role on small scales in biology and technology. For instance, biopolymers such as proteins or enzymes diffusing in biological cells intermittently bind to cellular membranes, or individual bacteria forming a biofilm on a surface use bulk excursions to efficiently relocate. Similarly the exchange between a liquid phase with a solid surface is an important phenomenon in the self-assembly of surface layer films and is a ubiquitous process in emulsions. This bulk-mediated surface diffusion, schematically shown in Fig. 1, was previously analyzed in terms of scaling arguments and simulations bychuk (); bychuk1 (); revelli (); fatkullin (), and was unveiled in field cycling NMR experiments in porous glasses stapf (). Moreover, effects of bulk-surface interchange were reported on proton transport across biological membranes membrane (). Recent studies are concerned with effects of bulk-surface exchange on reaction rates in interfacial systems benichou () and with surface diffusion of coppper atoms in nanowire fabrication nanowire ().

The remarkable finding of the bulk mediated surface diffusion model is that the effective surface motion is characterized by a Cauchy propagator bychuk (); bychuk1 (); revelli (); fatkullin ()

 ns(r,t)≃c1/2t2π(r2+ct2)3/2, (1)

where denotes a scaling property ignoring multiplicative constants, and is a dimensional factor. The associated stochastic transport is of superdiffusive nature bychuk (); bychuk1 (); revelli (); fatkullin (); stapf (); report (),

 ⟨r2(t)⟩s≃t3/2. (2)

Here we present a strictly analytical approach to this process. Our findings corroborate the previous scaling results for superdiffusion, however, we also derive the cutoffs to this behavior: at sufficiently long distances, the Cauchy propagator turns over to a Gaussian wing. Moreover at longer times the effective surface diffusion becomes subdiffusive, due to the fact that the particle spends less and less time on the surface. Normalized to the time-dependent surface coverage the effective surface diffusion turns over from superdiffusion to normal diffusion.

## Ii Derivation from a discrete random walk process

We start with a derivation of the coupling between surface and bulk in a discrete random walk process along the coordinate perpendicular to the surface (Fig. 1). Let with denote the number of particles at site of this one-dimensional lattice with spacing . The number of particles on the surface at lattice site are termed . The exchange of particles is possible only via nearest neighbor jumps. Each jump event between bulk sites is associated with the typical waiting time . For the exchange between the surface and site we then have the following law

 dN0(t)dt=12τN1−1τdesN0, (3)

where is the characteristic time for desorption from the surface. The probability of adsorption to the surface here is one, i.e., a particle adds to the surface population automatically when moving from site 1 to site 0. The exchange site at between the surface at and the next bulk site is governed by the balance relation

 dN1(t)dt=1τdesN0−1τN1+12τN2. (4)

Finally, the bulk sites are governed by equations of the form

 dN2(t)dt=12τN1+12τN3−1τN2, (5)

etc. Let us define the number of “bulk” particles at the surface site through

 N0≡2ττdesN0. (6)

This trick allows us to formulate the exchange equation also for site in a homogeneous form. Namely, from Eq. (3) we have

 dN0(t)dt=12τ(N1−N0). (7)

Moreover, from Eqs. (5) we find

 dNi(t)dt=12τ(Ni−1−2Ni+Ni+1), (8)

for .

Let us now take the continuum limit. For that purpose we make a transition from as the number of surface particles, and for the bulk concentration of particles. The factor may be viewed as the lattice constant measuring the distance between successive sites along the axis. As the bulk density is a function of while the surface particles are all assembled at , we need this length scale to match dimensionalities between and which are and . Expansion of the right hand side of Eq. (7) yields the surface-bulk coupling

 ∂ns(t)∂t=a2τa∂nb(z,t)∂z∣∣∣z=0. (9)

Similarly from Eq. (8) we obtain the bulk diffusion equation

 ∂nb(z,t)∂t=12τa2∂2nb∂z2. (10)

Finally, the boundary condition

 a2τnb(z,t)∣∣∣z=0=1τdesns (11)

stems from our definition (6). This completes the description of the particle exchange between surface and bulk, as well as the diffusion of bulk particles along the axis.

We now turn back to the full three-dimensional problem and add to above equations the directions (see Fig. 1). As the motion in the perpendicular direction is fully independent, we may simply adjust the Laplacian to three dimensions, ending up with the diffusion equation

 ∂nb(x,y,z,t)∂t=Db(∂2∂x2+∂2∂y2+∂2∂z2). (12)

with diffusivity . For the surface diffusion along the and coordinates with , we end up with the two-dimensional diffusion equation with diffusivity , plus the bulk-surface exchange term:

 ∂ns(x,y,t)∂t=Ds(∂2∂x2+∂2∂y2)ns+Db∂nb∂z∣∣∣z=0. (13)

These two equations are valid in the range and , and are supplemented by the initial condition

 ns(x,y,t)∣∣t=0=n0δ(x)δ(y) (14)

indicating that initially the particles are all concentrated on the surface at . Moreover we observe the boundary condition for ,

 nb(x,y,z,t)∣∣z=0=μns(x,y,t), (15)

and the bulk initial condition

 nb(x,y,z,t)∣∣t=0=0. (16)

In what follows, for simplicity of notation we use a unit initial concentration .

The surface and bulk diffusivities can be expressed in terms of the lattice constant and the typical waiting times between jumps in the bulk, , and on the surface, , through

 Db≡a26τ, (17)

and

 Ds≡a24τs, (18)

respectively. In the continuum limit, both and (or ) tend to zero such that the diffusion constants remain finite. In many physical systems bulk diffusion is considerably faster, such that . Above we also introduced the coupling parameter

of physical dimension . Small values of at fixed and correspond to slow bulk-surface exchange.

For consistency we derive the overall number of particles. To this end we define the number particles on the surface, , and in the bulk, through the relations

 Ns(t)=∫∞−∞dx∫∞−∞dyns(x,y,t) (20)

and

 Nb(t)=∫∞0dz∫∞−∞dx∫∞−∞dyns(x,y,t). (21)

From integration of Eq. (13) we find

 dNs(t)dt=Db∂∂z∫∞−∞dx∫∞−∞dynb(x,y,z,t)∣∣∣z=0. (22)

Similarly, Eq. (12) yields

 dNb(t)dt = Db∫∞0dz∂2∂z2∫∞−∞dx∫∞−∞dynb(x,y,z,t) (23) = −Db∂∂z∫∞−∞dx∫∞−∞dynb(x,y,z,t)∣∣∣z=0.

Combination of Eqs. (22) and (23) produces

 ddt(Ns(t)+Nb(t))=0. (24)

Thus the overall number of particles is conserved, as it should be.

## Iii Solution of the coupled diffusion problem

To solve the set of coupled equations (13) and (12) for the specified initial and boundary value problem, we start by defining the two-dimensional Green’s function

 Gs(x,x′,y,y′,t)≡Gs(x−x′,y−y′,t)=1√4πDstexp(−(x−x′)2+(y−y′)24Dst). (25)

Then, the solution of Eq. (13) becomes

 ns(x,y,t) = ∫dx′∫dy′Gs(x,x′,y,y′,t)ns(x′,y′,t=0) (26) +∫t0dt′∫dx′∫dy′Gs(x,x′,y,y′,t−t′)(Db∂nb(x′,y′,z,t′)∂z)z=0.

With initial condition (14) we thus find

 (27)

where

 Gs(kx,ky,s)=1s+Db[k2x+k2y] (28)

is the Fourier-Laplace transform of the surface Green’s function . Here and in the following we denote the Laplace and Fourier transform of a function by explicit dependence on the image variables, that is,

 f(s)=L{f(t);t→s}=∫∞0f(t)e−stdtandg(k)=F{g(x);x→k}=∫∞−∞g(x)eikxdx. (29)

The bulk particle density according to Eq. (12) is given by the formal expression jaeger ()

 nb(x,y,z,t)=z(4πDb)3/2∫t0dt′∫dx′∫dy′μns(x′,y′,t′)(t−t′)5/2exp(−(x−x′)2+(y−y′)2+z24Db(t−t′)). (30)

From this expression one can indeed show that, despite the factor , the coupling equation (15) is fulfilled, see Appendix B. Now, Eq. (27) requires the derivative of expression (30) with respect to , evaluated at . To find that expression, we first differentiate

 ∂nb∂z=μ(4πDb)3/2∫t0dt′∫dx′∫dy′ns(x′,y′,t′)(t−t′)5/2exp(−(x−x′)2+(y−y′)2+z24Db(t−t′))[1−z22Db(t−t′)], (31)

and then calculate its Fourier-Laplace transform

 L{F{Db(∂nb∂z)z=0}}=μDb√4πDbns(kx,ky,s)L{exp(−Db[k2x+k2y]t)t3/2}z=0, (32)

due to the convolution nature of expression (31). The Laplace transform is evaluated by help of the shift theorem, yielding

 L{exp(−Db[k2x+k2y]t)t3/2}=L{1t3/2}s→s+Db[k2x+k2y]=−√4π(s+Db[k2x+k2y]). (33)

We thus finally obtain

 L{F{Db(∂nb∂z)z=0}}=−μD1/2b√s+Db(k2x+k2y)ns(kx,ky,s). (34)

Insertion of this equation into expression (27) delivers the solution for the surface density in Fourier-Laplace space,

 ns(kx,ky,s)=Gs(kx,ky,s)1+μD1/2bGs(kx,ky,s)√s+Db[k2x+k2y], (35)

where the Fourier-Laplace transform of the Green’s function was defined in Eq. (28). After some transformations we arrive at the exact closed form expression

 ns(kx,ky,s)=1s+Ds[k2x+k2y]+χ√s+Db[k2x+k2y], (36)

with the rescaled coupling parameter

of dimension . Relation (36) is the main result of this work, and we now consider the consequences to the surface motion effected by the bulk mediation.

## Iv Effective surface behavior

We first determine the number of particles on the surface,

 Ns(t)=∫∞−∞∫∞−∞ns(x,y,t)dxdy, (38)

whose Laplace transform is

 Ns(s)=n(kx,ky,s)∣∣kx=ky=0=1s+χs1/2. (39)

Inverse Laplace transformation then yields the exact expression

 Ns(t)=eχ2terfc(χt1/2), (40)

where we use the complementary error function

 erfc(z)=2√π∫∞ze−ξ2dξ=1−erf(z). (41)

At short times , this leads to the initial decay

 Ns(t)∼1−2χ√πt1/2 (42)

of the number of surface particles, eventually turning into the long time behavior

 Ns(t)∼1χ√πt. (43)

The asymptotic decay stems from the returning dynamics to the origin of a Brownian motion along the coordinate, i.e., it is proportional to the normalization factor of a one-dimensional Brownian motion. The additional prefactor rescales time with respect to the efficiency of the surface-bulk exchange.

We now turn to the surface dynamics, as quantified by the effective mean squared displacement along the surface. In the Laplace domain,

 ⟨r2(s)⟩s = −∇2kx,kyn(kx,ky,s)∣∣kx=ky=0 (44) = −[∂2∂k2+1k∂∂k]n(kx,ky,s)∣∣kx=ky=0 = 4Dss(√s+χ)2+2χDbs3/2(s1/2+χ)2,

where . From this expression we obtain the limiting behaviors at short and long times. Thus, we observe that the short time limit in Laplace domain corresponds to , such that

 ⟨r2(s)⟩s∼4Dss2+2χDbs5/2. (45)

This translates into the asymptotic time evolution

 ⟨r2(t)⟩s∼4Dst(1+23√πDbDs[tχ2]1/2). (46)

As long as the ratio is sufficiently large, there is a superdiffusive component winning over the normal surface diffusion proportional to . This is exactly the famed bulk mediated superdiffusion originally obtained from scaling arguments by Bychuk and O’Shaugnessy bychuk (); bychuk1 (). Note that this diffusional enhancement is accompanied by an almost constant number of surface particles, compare Eq. (42).

Conversely, at long times (or in Laplace domain) we find

 ⟨r2(s)⟩s∼4Dsχ2s+2Dbχs3/2, (47)

corresponding to the temporal behavior

 ⟨r2(t)⟩s∼4Dsχ2(1+1√πDbDs[tχ2]1/2). (48)

Somewhat surprisingly, at sufficiently long times the bulk contribution to the effective mean squared displacement dominates over the surface contribution for arbitrary ratio , giving rise to subdiffusive behavior. This subdiffusion occurs due to the ongoing loss of surface particles into the bulk, see Eq. (43).

Instead of considering the surface mean squared displacement we introduce the normalized effective surface mean squared displacement

 ⟨r2(t)⟩norms≡1Ns(t)⟨r2(t)⟩s. (49)

This quantity can be interpreted as the surface mean squared displacement covered by an individual particle that effectively stays on the surface and does not fully escape to the bulk. At long times this quantity has the limiting form

 ⟨r2(t)⟩norms∼4Dbt. (50)

The long time diffusion corrected for the number of escaping particles displays normal diffusion, albeit with the bulk diffusivity.

In the following we neglect contributions from the surface diffusion proportional to , in order not to overburden the presentation. The interesting behavior is due to the bulk mediation with weight . We quantify the motion in terms of fractional order moments, before embarking for the surface propagator.

### Fractional order moments

We now derive an exact expression for the -th order moments ()

 ⟨|r|q(t)⟩s=∫|r|qns(r,t)d2r. (51)

To this end we utilize the following integral on the plane:

 ∫(1−cos(k⋅r))dk|k|2+q=2π∫∞0(1−J0(kr))dkk1+q=2πrq∫∞0(1−J0(z))dzz1+q, (52)

where we used polar coordinates and corresponding to the two-dimensional vectors and . Thus we identify

 rq=K(q)∫(1−cos(k⋅r))dkk2+q, (53)

with the definition

 K(q)=(2π∫∞0(1−J0(z))dzz1+q)−1=2qπ2sin(πq2)[Γ(1+q2)]2. (54)

With this trick we can rephrase the th order moment (51) as follows,

 (55)

The integral over of the surface density is but the number of surface particles, such that

 ⟨rq(t)⟩s=K(q)∫d2kk2+q[Ns(t)−Re{ns(k,t)}], (56)

where we replaced the Fourier cosine transform of by the real part of the exponential Fourier transform.

With the Fourier-Laplace transform (36) of the surface propagator with set to zero, for the Laplace transform of the th order moment we obtain

 ⟨rq(s)⟩s = K(q)∫[ns(k=0,s)−ns(k,s)]dkk2+q=2πK(q)[1s+χ√s−1s+χ√s+Dbk2]dkk1+q (57) = 2πK(q)χs+χ√s∫∞0√s+Dbk2−√ss+χ√s+Dbk2dkk1+q = 2πK(q)χDbs+χ√s∫∞01(√s+√s+Dbk2)(s+χ√s+Dbk2)dkkq−1 = 2πK(q)χ2s+χ√s(Dbs)q/2∫∞01(χ+χ√1+y2)(√s+χ√1+y2)dyyq−1 = 2πK(q)χDq/2bs(1+q)/2(s−χ2){I1(q)−I2(q)},

where on the way we introduced the substitution . The two integrals are defined by

 I1(q)=∫∞0y1−q1+√1+y2dy (58)

and

 I2(q)=∫∞0y1−q√s/χ+√1+y2dy. (59)

Now we analyze the temporal behavior of the th order moment at short and long times.

#### iv.0.1 Long time behavior

At long times (or ) we see that

 I1(q)−I2(q)∼∫∞0(11+√1+y2−1√1+y2)dyyq−1<0, (60)

and

 ⟨rq(s)⟩s≃1s(1+q)/2, (61)

corresponding to the time evolution

 ⟨rq(t)⟩s≃t(q−1)/2. (62)

Taking normalization by the number of surface particles into account, we find

 ⟨rq(s)⟩norms=⟨rq(s)⟩sNs(t)≃tq/2. (63)

That is, at long times the surface diffusion exhibits normal scaling behavior.

#### iv.0.2 Short time behavior

The more interesting case is the short time behavior corresponding to the limit (or ). Here we consider three separate cases:

(i) The case : Since both integrals and converge, we may simply neglect . Then

 ⟨rq(s)⟩s≃s−(3+q)/2, (64)

such that after Laplace inversion we find

 ⟨rq(t)⟩s≃t(q+1)/2. (65)

(ii) The case : Now we should take into account both integrals,

 I1(q)−I2(q) ∼ (66) √sχ∫∞0y1−qdy(1+√1+y2)(√s/χ+√1+y2).

To estimate the main contribution from this difference, we split it into three parts, namely

 I1(q)−I2(q) ∼ (67) √sχ{∫10…dy+∫√s/χ1…dy+∫∞√s/χ…dy},

and evaluate each contribution separately. We find

 ∫10…dy≃χ√s, (68a) and then ∫√s/χ1…dy∼∫√s/χ11y√s/χdyyq−1∼(χ√s)q>χs. (68b) Finally, ∫∞√s/χdyy1+q≃(χ√s)q. (68c)

Thus the main contribution come from Eqs. (68b) and (68c), and

 I1(q)−I2(q)≃s(1−q)/2. (69)

With this estimate the Laplace transform of the th order moment, Eq. (57), has the leading order behavior

 ⟨rq(s)⟩s≃s(1−q)/2s(3+q)/2≃1s1+q. (70)

This corresponds to

 ⟨rq(t)⟩s≃tq (71)

in the time domain.

(iii) The case : This special case requires some care. We start with the substitution in Eqs. (58) and (59). Then the integrals become

 I1(1)=∫π/20dϕcosϕ(1+cosϕ), (72)

and

 I2(1)=∫π/20dϕcosϕ(1+√scosϕ/χ). (73)

We can now rewrite Eq. (57) in the form

 ⟨r(s)⟩s = 2πK(1)χD1/2bs(s−χ2) (74) × (√sχ∫π/20dϕ1+√scosϕ/χ−1),

where we used prudnikov ()

 ∫π/20dϕ1+cosϕ=1. (75)

For the integral in the parenthesis of Eq. (74) becomes prudnikov ()

 ∫π/20dϕ1+√scosϕ/χ=2√1−s/χ2arctan√1−√s/χ1+√s/χ. (76)

Thus, for Eq. (74) we find

 ⟨r(s)⟩s = 2πK(1)χD1/2bs(s−χ2) (77) ×(2√sχ1√1−s/χ2arctan√1−√s/χ1+√s/χ−1),

so that we find the behavior

 ⟨r(s)⟩s∼D1/2bχs. (78)

After Laplace inversion,

 ⟨r(t)⟩s∼D1/2bχ,and⟨r(t)⟩norms∼(Dbt)1/2 (79)

at long times, , consistent with Eq. (63). Conversely, for we employ prudnikov ()

 ∫π/20dϕ1+√scosϕ/χ = 1√s/χ2−1 (80) ×ln√s/χ2−1+√s/χ−1√s/χ2−1+1−√s/χ.

Thus we find

 (81)

We thus obtain the limiting form at large

 ⟨r(s)⟩s∼χD1/2bs2ln(2√sχ). (82)

Back-transformed this results in

 ⟨r(t)⟩s≈⟨rq(s)⟩norms∼tlnt (83)

at short times .

Summarizing our results for th order moments, at short times the effective surface diffusion exhibits anomalous scaling: the th order moment of the radius scales like for and for , while the first moment includes a logarithmic contribution, , consistent with the earlier results in Ref. bychuk1 (). At long times the th order moments scale normally with time, proportional to .

## V Surface propagator

We now turn to the behavior of the surface propagator. If we neglect surface diffusion (i.e., ) the surface propagator from Eq. (36) becomes

 ns(k,s)=1s+χ√s+Dbk2. (84)

We perform an inverse Laplace transformation along the Bromwich path:

 ns(k,t) = ∫Brexp(st)s+χ√s+Dbk2ds2πi=∫Brexp(sχ2t)s+√s+κ2ds2πi=e−κ2χ2t∫Brexp(sχ2t)s+√s−κ2ds2πi (85) =e−κ2χ2t∫Brexp(sχ2t)(√s+12+√κ2+14)(√s+12−√κ2+14)ds2πi =exp(−κ2χ2t)2√κ2+1/4⎛⎜ ⎜⎝∫Brexp(sχ2t)√s+12−√κ2+14ds2πi−∫Brexp(sχ2t)√s+12+√κ2+14ds2πi⎞⎟ ⎟⎠ =exp(−κ2χ2t)2√κ2+1/4⎛⎜⎝L−1{1√s−b;s→t′}b=√κ2+1/4−1/2,t′=tχ2−L−1{1√s+a;s→t′}a=1/2+√κ2+1/4,t′=tχ2⎞⎟⎠

where on the way we introduced the substitutions and . With

 L−1{1√s+a;s→t}=1√πt−aea2terfc(a√t) (86)

we finally arrive at the Fourier transform of the surface propagator,

 ns(k,t) = 12√κ2+1/4{[12+√κ2+14]exp([12+√κ2+14]tχ2)erfc([12+√κ2+14]χ√t) (87) +(√κ2+14−12)exp([12−√κ2+14]tχ2)