Bulk-mediated surface diffusion: non-Markovian desorption dynamics

Here we analyse the dynamics of adsorbed molecules within the bulk-mediated surface diffusion framework, when the particle's desorption mechanism is characterized by a non-Markovian process, while the particle's adsorption as well as its motion in the bulk is governed by Markovian dynamics. We study the diffusion of particles in both semi-infinite and finite cubic lattices, analysing the conditional probability to find the system on the reference absorptive plane as well as the surface dispersion as functions of time. The results are compared with known Markovian cases showing the differences that can be exploited to distinguish between Markovian and non-Markovian desorption mechanisms in experimental situations.


Introduction
Since the pioneering works of Einstein [1], Smoluchovski [2] and Langevin [3], the study of diffusion processes and their application has pervaded all areas of science and technology. In addition to some well-known encyclopaedic works [4,5], over the last few decades, the literature on diffusion processes has covered from some formal mathematical aspects [6], through diffusion-limited reactions [7,8] and anomalous diffusion or diffusion in disordered systems [9], to heat conductivity in solids [10] as well as other diffusion problems in material science (ranging from polymers to hydrogen in metals [11,12]). Other problems also studied within the framework of diffusion processes include combustion [13], applications in cardiology [14,15] and biology [16]- [19], in self-organization and pattern-forming systems [20]- [22], up to applications in economic and sociological contexts [23]- [25].
Among the many problems studied in material science, the dynamics of adsorbed molecules at an adsorbing surface is a fundamental issue in interface science [26]- [31] and is crucial to a number of emerging technologies [26,32]. Its role is central to phenomena as diverse as foam relaxation [33] and the evolution of blood protein deposits [34]. Recently, the mechanism called bulk-mediated surface diffusion has been identified and explored. The importance of bulksurface exchange in relaxing homogeneous surface density perturbations is experimentally well established [35]- [39]. This mechanism typically arises at interfaces separating a liquid bulk phase and a second phase which may be either solid, liquid or gaseous. Whenever the adsorbed species are soluble in the liquid bulk, adsorption-desorption processes occur continuously. These processes generate surface displacement because the desorbed molecules undergo Fickian diffusion in the liquid's bulk, and are later re-adsorbed elsewhere. When this process is repeated many times, an effective diffusion results for the molecules on the surface.Adsorption at the solidliquid interfaces arises, for instance, in many biological contexts involving protein deposition [40]- [42], in solutions or melts of synthetic macromolecules [43]- [46], in colloidal dispersions [47] and in the manufacture of self-assembly mono-and multi-layers [26,32,48,49].
Usually, the studies performed in this type of systems are carried out within the framework of a Master Equation scheme [50]- [52], where the particle's motion through the bulk and the adsorption-desorption processes are Markovian. In two recent papers [53,54], we have shown some of the most important features of this phenomenon. In particular, by studying the variance r 2 (t) plane of the position r ≡ (x, y) on the interface, or the conditional probability P(z = 1, t) ≡ x,y P(x, y, z = 1; t|0, 0, 1; t = 0) to find the particle on the interface at time t, if it was initially at (0, 0, 1), through analytical and numerical approaches, the following results were obtained. For the case of a semi-infinite or unbounded bulk [53] • for t → ∞, the effective diffusion on the interface (first layer of the lattice) is always subdiffusive (the variance of the position grows as t 1/2 ) regardless of the desorption rate δ. Similarly, the probability to find the particle on the interface at time t decays as t −1/2 , independently of δ; • at finite times, the growth of the variance can be fitted by a t law. The exponent depends on the range of time considered and the values of the adsorption and diffusion constants, increasing rapidly as δ decreases and saturating at a value compatible with the one reported in [50,51]; • an effective continuous-time random walk (CTRW) description (without conservation of probability) was derived on the interface.
For the finite or bounded bulk case [54], we have investigated the transition from multilayered to unbounded bulk regime, and found that • there exists an optimal number of layers that maximizes r 2 (t) plane on the interface (which is a measure of the effective diffusivity); • and up to about that thickness, the long-time effective diffusivity on the interface has normal character, and crosses over abruptly towards a sub-diffusive behaviour as the number of layers increases further.
It is worth remarking here that for an arbitrary (finite) number of layers, the Laplace transform usually cannot be analytically inverted. This forced us to apply numerical inversion methods whose efficacy has been tested-with excellent results-not only against analytically solvable cases, like the bilayer one, but also against Monte Carlo simulations.
Here we address the dynamics of adsorbed molecules when the particle's desorption mechanism is characterized by a non-Markovian process, while the particle's adsorption and its motion in the bulk are governed by a Markovian dynamics. A non-Markovian desorption process can occur when the surfaces contain 'deep traps', capture and re-emission from surfaces that contain sites with several internal states such as the 'ladder trapping models', proteins with active sites deep inside its matrix, etc [55,56].
When the non-Markovian process is present it is necessary to resort to Generalized Master Equations. These equations are characterized by a 'memory kernel'and may be related univocally with a CTRW scheme (for instance, see the paper by Haus and Kehr [55]).
The main goal of this work is to study the influence of the non-Markovian desorption dynamics on the effective diffusion process at the interface z = 1. For that purpose, we calculate the temporal evolutions of the variance ( r 2 (t) plane ) and P(x, y, z = 1; t|0, 0, 1, t = 0), which (as indicated above) is the conditional probability of finding the particle on the surface at time t since the particle arrived there at t = 0 which, in what follows, we indicate by P(z = 1, t).
In the next section we formally present the model, in terms of a Generalized Master Equation which describes the particle's dynamics through the bulk and its desorption. In order to make the paper self-contained, and for future comparison, we start with a brief review of the Markovian desorption case, and then indicating the differences that are introduced in the formalism by the non-Markovian desorption process. Next, we present some analytical and numerical results for the non-Markovian desorption case in both situations: for infinite as well as for finite systems. In the following section, and for an infinite system, we discuss the behaviour when the waiting time density has a 'long time tail'. Finally, in the last section we discuss the results and present some conclusions.

Markovian case
Here we start with a brief review of the model we used, as applied in [53] to the Markovian and infinite cases. We consider the problem of a particle making a random walk in a semi-infinite cubic lattice. The 'bulk' is unbounded in the z direction, where the particles can move from z = 1 to z = ∞ in the infinite case (from z = 1 to z = L, with L the number of layers in the finite or bounded case), while the x and y directions remain unbounded. The position of the walker is defined by a vector r whose components are denoted by the integer numbers n, m, l corresponding to the directions x, y and z, respectively.
The probability that the walker is in (n, m, l) at time t, given that it was at (0, 0, l 0 ) at t 0 , which we indicate by P(x, y, z = 1; t|0, 0, 1, t 0 = 0) ≡ P(n, m, l; t), satisfies the following set of coupled master equationṡ where α, β and γ are the transition probabilities per unit time in the x, y and z directions, respectively, and δ is the desorption probability per unit time from the boundary plane defined by z = 1.
Taking the Fourier transform with respect to the x and y variables and the Laplace transform with respect to the time t in the above equations, we obtain Above, we have used the following definitions: where L indicates the Laplace transform of the quantity within the brackets, and It is possible to write equation (2) in matrix form as where the square matrixG has components In equation (5),Ĩ is the identity matrix andH a three-diagonal matrix with the following form: and C is defined as In order to find the solution of equation (5), in [53] we have decomposed theH matrix into different contributions and exploited a Dyson-like approach in order to obtain analytical expressions (in the Laplace domain) for • r 2 (t) plane , the variance of the probability distribution at time t on the plane z = 1, • P(z = 1, t), the probability that the particle is in the plane z = 1 at time t.
The variance r 2 (t) plane , is a direct measurable experimental magnitude [50,51], which measures the spreading of particles over the z = 1 plane. Once P(m, n, l = 1; t) = P(m, n, l = 1; t|0, 0, l 0 = 1; t = 0) is known, the variance is calculated as where the symmetry properties for the diffusion along the x and y axes, that is x(t) = y(t) = 0, were used. As indicated in [53], it is important to remark that the conservation of particles in the plane is not satisfied. For P(z = 1, t), it was shown [53] that the Laplace transform is In [53], we have used a numerical procedure to obtain the inverse Laplace transform of both r 2 (t) plane and P(z = 1, t), and compared these theoretical results with Monte Carlo simulations and found an excellent agreement.
Clearly, for finite systems, that is for a finite number of layers, the set of equations indicated by equations (1) and (2), reduces to a finite number (up to L, the number of layers) of coupled equations. We will not repeat it here but instead refer to [54].

Non-Markovian case
We turn now to the non-Markovian desorption. We use the same framework as before. Now P(n, m, l; t) satisfies the following Generalized Master Equatioṅ where α, β and γ, as before, are the transition probabilities per unit time through the bulk in the x, y and z directions, respectively, and K(t) represents the memory kernel for desorption at all sites over the z = 1 surface (that is for (n, m, l = 1)).
Taking again the Fourier transform with respect to the x and y variables and the Laplace transform with respect to the time t in the above equations, we obtain Here we have used the same definitions as in equations (3) and (4). Clearly, the above equations for G(k x , k y , l; s)] are similar to equation (2), with K(s) instead of δ and, therefore, all results obtained in [53] remain valid for a non-Markovian dynamics when δ is replaced by K(s). Namely, the Laplace transform of the variance r 2 (s) plane = L[ r 2 (t) plane ] can be found as [53,54] By using equations (5) and (12), r 2 (t) plane turns out to be the ratio of two functions of s where and For P(z = 1, t), according to the result shown in [53], we have We reiterate that, as indicated in [53], the conservation of particles on the plane is not satisfied.
To conclude this section, we address the relation between ψ(t), the waiting time density, as defined in a CTRW scheme, and the memory kernel of equation (10) (see the paper by Haus and Kehr [55], etc). In the Laplace domain, this relation reads

Analytical results and Monte Carlo simulations
In this section we show the results obtained analytically and by Monte Carlo simulations. In all simulations we have fixed the following parameters: α = β = γ = 1. The simulation results were obtained by averaging over 2 × 10 6 realizations.
In what follows, and in order to describe the desorption dynamics at the surface, we have used a family of waiting time densities that, since it was first introduced by Scher and Lax [57] to describe the frequency dependence of the electric conductivity in disordered solids when transport is due to impurity hopping, has been extensively exploited in modelling non-Markovian situations (see also [58]). The reason for its wide use are: the versatile functional form and its simplicity that allows to take into account a spread of transition rates in a controllable way [57]. When only one transition rate is present a Markovian description is reobtained: the memory kernel becomes a Dirac δ-function; when the spread is very wide, steps occur at quasi-periodic intervals of time (see below). The adopted function, which has a finite first moment, is where a is a positive integer and (a) the gamma or factorial function. It is worth remarking here two important facts about this family of functions. Firstly, as can be seen from equation (17) corresponds to the Markovian case; a = 1 to the non-Markovian case, while the parameter θ is the 'average desorption rate'. Secondly, the mean value of these waiting time densities is Such a result means that the mean-time value of this family of functions does not depend on the parameter a, but is the only function of the desorption rate. In figure 1 we depict some typical curves for this family of functions. Figure 2 depicts the temporal evolution for the conditional probability P(z = 1, t), for the Markovian and non-Markovian cases (a = 75 and θ = 0.01) of a system whose bulk is unbounded. From this figure we can observe two temporal regions: a transient region which ranges from t = 0 to 1000 and an asymptotic region where the behaviour approaches that of  the Markovian case. An important feature to be detached here is the arising of the damped oscillations in the transient region for the non-Markovian case. This oscillatory behaviour is due to the waiting time density functions used, as will be explained later.

Infinite systems
In figure 3 we depict the results for the variance r 2 (t) plane , for the same situation (a = 75 and θ = 0.01) and for the same Markovian and non-Markovian cases. We can observe two important features. The first one is that, in the non-Markovian regime, the system shows a delay in the beginning of the spreading. The second feature is that the system presents an oscillatory-like behaviour, with the oscillations attenuating as the time increases. This fact is still more apparent in the inset of the figure, where we show the temporal evolution of the spreading velocity (V Spread ).
It is worth remarking here that oscillations appear only in the non-Markovian case and are due to a particular behaviour of the family of waiting time densities defined in equation (17). When the markovianicity parameter tends to infinity, a → ∞, it is well known that ψ(t) tends towards a Dirac δ function: ψ(t) → δ(t − θ −1 ). This fact implies a kind of 'periodicity', as there is a desorption after each elapsed t θ −1 period. This is the origin of the oscillatory behaviour in both, P(z = 1, t) and r 2 (t) plane . It is also the origin of the initial delay. The figures also show an excellent agreement between the theoretical and Monte Carlo simulation results.

Finite systems
So far we have treated the infinite problem, observing the different behaviour for the Markovian and non-Markovian desorption processes, and the effects that they produce in the temporal evolution of the P(z = 1, t) and r 2 (t) plane . Now we will consider the finite systems, that is systems that are bounded in the normal direction to the observation surface, while the other directions (parallel to the plane) remain unbounded. We will show some new features and confirm others for this physical phenomenon. The temporal evolutions of r 2 (t) plane and P(z = 1, t) are depicted in figure 4 for the case of a bilayer system. The parameters used in this figure are similar to those used in figure 3.
Here we can observe that the oscillations arising in the evolution are larger than those produced in the infinite system (see both figures 2 and 3). This fact is corroborated in the inset of the figure. There, we have depicted the evolution of P(z = 1, t), observing that the non-Markovian system converges to the Markovian system but now to the values corresponding to the finite systems. The apparent differences between the oscillations on the r 2 (t) plane evolutions for the finite and infinite bulk cases are better depicted in figure 5, where we have shown the V Spread for three finite systems with different number of layers. As can be seen, a larger number of layers produce a smaller amplitude of the oscillation.
Another important feature is shown in figure 6, where we have depicted the evolution of P(z = 1, t) for the bilayer case and for two different values of a (a = 50 and 700). From the figure we can see that the effect produced by the non-Markovian desorption behaviour is a change on the oscillation's amplitude. The inset of the figure shows the dependence of the amplitude of the first oscillation peak (which we indicate by ) as a function of a, the markovianicity parameter. It is important to remark the way in which the non-Markovian behaviour affects the amplitude of oscillation.
From this figure, there is another point to be remarked upon: it is that ω, the oscillation frequency, remains unperturbed. Moreover, the non-Markovian behaviour does not affect the frequency of the oscillation. This fact is again due to the particular behaviour of the family of waiting time densities defined in equation (17). This frequency directly depends on the mean residence time (θ −1 ) on the surface. A simple analysis shows that a linear relation ω ≈ θ is fulfilled. This fact is clearly seen in figure 7. Finally, in figure 8 we depict P(z = 1, t) as a function of t for the bilayer case when the adsorption rate γ is varied. As γ increases, the oscillation amplitude decreases, and in addition the particles remain more time adsorbed on the surface. The inset shows again the linear relation between ω and θ, which results to be independent of γ.

Behaviour for 'long tail' waiting time densities for infinite systems
Here we analyse the asymptotic-large t-behaviour of r 2 (t) plane and P(z = 1, t). As usual, in order to obtain such a behaviour we analyse the s 1 limit, for instance from equations (13)- (15). We assume that ψ(t), the waiting time density for desorption, has an infinite first moment. This means that, when s 1, ψ(s) ∼ 1 − Bs ν with 0 < ν < 1. Consequently, in this limit K(s) ∼ (1/B)s 1−ν . From equations (13)-(16), we obtain and The behaviour of the denominators in equations (19) and (20) depends on the values adopted by the exponent ν. In the different ranges we have the following behaviour for r 2 (t) plane and From these expressions, and by means of Tauberian theorems [58], the behaviour for large t results and It is worth remarking here that equations (23) and (24) represent the behaviour of r 2 (t) plane and P(z = 1, t), in the long time limit, for any non-Markovian desorption waiting time density function with a long tail in time.
From equation (23) it can be deduced that the asymptotic sub-diffusive regime arises for all values of ν (0 < ν < 1) except for the case ν = 1 2 , where normal diffusion take place. This occurs due to the 'competition' between adsorption and Fickian diffusion across the bulk. We have shown in [53] that the probability density to return for the first time to the plane z = 1 is proportional to t −3/2 for large times; it is the behaviour of the desorption waiting time density function with ν = 1 2 that produces the particular value of this parameter. The ladder trapping model [55] with infinite number of internal states is an important example, where the desorption waiting time density is characterized by a parameter ν = 1 2 . In figure 9 we show the exponent of t, r 2 (t) plane ∼ t , corresponding to equations (23). Although for infinite systems 1 2 1 we have recently shown [60] that for finite-or infinitebiased systems this exponent may assume values between 0 and 1. The following figure shows the behaviour of P(z = 1, t) and r 2 (t) plane for the following desorption waiting time density We remark that we have fixed the parameters α, β and γ equal to one. We show in figures 10 and 11 the behaviour of P(z = 1, t) and r 2 (t) plane , respectively, as functions of t, for φ = 0.1 and different values of ν = 0.1. These figures clearly depict the asymptotic behaviour indicated above.

Conclusions
We have studied here the evolution of particles making an effective diffusion on a surface. The diffusion is actually performed across the bulk surrounding the surface, resulting in the socalled bulk-mediated surface diffusion phenomenon. Usually, the models proposed are based on a Markovian process. The main feature of this work was to present an analytical model for non-Markovian desorption from the surface. For the bulk that surrounds the surface we have considered the two possibilities, it being semi-infinite or finite, and the particles undergoing a Markovian motion there. It is worth remarking here that, even though we have only considered diffusion on the bulk but not on the surface, the model also allows us to treat the general case with simultaneous diffusion in bulk and surface.
We observed the influence of the non-Markovian desorption dynamics on the effective diffusion process at the interface z = 1 by calculating analytically, in the Laplace domain, the temporal evolutions of the variance ( r 2 (t) plane ) and the conditional probability of being on the surface at time t since the particle arrived there at t = 0 (P(z = 1, t)). We have chosen a family of non-Markovian desorption waiting time densities which has a finite first moment and compared the analytical results for r 2 (t) plane and P(z = 1, t) with Monte Carlo simulations obtaining an excellent agreement for both: semi-infinite and finite systems. For these particular desorption waiting time densities, we can establish two regions based on the ratio θ/γ, calling them strong adsorption (when θ/γ 1) and weak adsorption (when θ/γ 1). In the strong adsorption limit, we can observe a transient regime and an asymptotic one. The main feature in the transient regime is that the temporal evolution is characterized by damped oscillations whose frequencies are in direct relation to θ, the average desorption rate. It is worth remarking here that the frequency does not depend on the markovianicity degree, but the effect of this parameter appears on the amplitude of the signal. In addition, when we consider a system with a finite bulk (that is, with a finite number of layers), the oscillation increases its amplitude. In the asymptotic regime, the conditional probability of the non-Markovian system tends to the Markovian system. Finally, the oscillatory effects disappear in the weak adsorption region for both semi-infinite and finite systems, where the behaviour approaches the Markovian one.
We have also analytically studied the effective diffusion processes when the waiting time density for desorption has an infinite first moment. We have deduced that asymptotic sub-diffusive regime appears for all values of the parameter ν, which characterize the behaviour of the 'tail' of ψ(t) (0 < ν < 1), except for the case ν = 1 2 where normal diffusion takes place. When finite-or infinite-biased systems are considered [60], the dispersion increases and gives rise to a slower growth of the variance with time compared with that of the infinite unbiased systems. Such lower values of the exponent are in correspondence with some experimental results [61,62].
Finally, it is worth remarking here an important aspect of the present approach. Through the above results we have shown that the behaviour of r 2 (t) plane and of P(z = 1, t) are strongly dependent on the desorption mechanism. As the effective dispersion and the percentage of particles that remain in z = 1 are of measurable magnitudes [51] they may be used to investigate the characteristic and fundamental parameters of the desorption processes.