Kinetics and energetics of reflected diffusion process with time-dependent and space-homogeneous force

We present the exact analysis of a spatially restricted one-dimensional diffusion process with a time-dependent, spatially linear potential and a reflecting boundary. Using matching conditions together with the known solution for the unbounded diffusion in the time modulated potential, we have derived a new integral equation, whose solution yields the Green function for the restricted diffusion. Applying the general scheme, we give the numerical analysis of the diffusion in a symmetrically oscillating force field superimposed on the time-independent component. The latter component alone guarantees the approach to Gibbs equilibrium with exponential probability density. We calculate both the kinetic and the energetic characteristics of the emerging non-equilibrium isothermal process and discuss their dependence on the model parameters.


Introduction
Hypothetically, assume that the intensity of the uniform gravitational field acting on the isothermal atmosphere above the Earth surface is periodically modulated. How would the air pressure at sea level change? How does the centre of gravity of the air column move, i.e. how does the potential energy of the air change? How much work must be done to keep the system in the stationary non-equilibrium regime? These fairly oversimplified questions guide us to an interesting problem within diffusion theory. Its solution forms the main objective of this paper.
The dynamical behaviour of an overdamped Brownian particle acted upon by the thermal force and moving in a fixed potential landscape is a well understood classical problem [1]- [3]. However, in the last decade, the diffusion dynamics has been re-examined in systems in which the potential depends on time, the modulation being due to additional (deterministic [4], stochastic [5], or internal self-organized [6]) dynamical mechanisms [7]. The new achievements have substantially broadened the field (cf the two recent Focus Issues [8,9]).
The diffusion dynamics in time-dependent potentials plays a central role in the phenomenon of stochastic resonance [10], in physics of Brownian motors [4,11,12] and in the recent discussion concerning the energetics of the diffusion process [13]-these recent reviews discuss history, applications and existing literature in the domain. Diffusion in a time-dependent potential can be regarded as an example of an isothermal irreversible process. Identifying the work done on the system by the external agent and the heat exchange with the heat bath [14,15] one immediately enters the discussion of the famous Clausius inequality between the irreversible work and the free energy. From this perspective, the work done on the system has been associated with individual realizations (trajectories) of the diffusive motion, i.e. the work itself is treated as a random variable whose mean value enters the thermodynamical considerations. An important achievement in the field is the discovery of the new fluctuation theorems, which generalize the Clausius identity in giving the exact mean value of the exponential of the work. This Jarzynski identity [16,17] enables one to specify of the free energy difference between two equilibrium states. This is done by repeating real time (i.e. non-equilibrium) experiment and measuring the work done during the process. The identity has been recently experimentally tested [18].
At the same time, there is only a limited number of non-trivial exactly solvable diffusion models with a time-dependent potential, cf for example [19]- [22]. This is an unsatisfactory situation because an analytical insight into the solution represents a considerable advantage. Indeed, the above mentioned physical effects are mostly conditioned by the comparability of several time-and space-scales. Their theoretical description precludes the usual perturbational approaches.
Our main interest here concerns a physically relevant situation with a time-dependent potential which, up to our knowledge, does not yet seem to have been investigated in the literature. The central idea is to use the known exact solution for the time-dependent potential and incorporate in it the reflecting boundary. The spatial restriction of the original diffusion problem implies a new physically interesting regime. The flexibility of the method itself rests in the fact that the original spatially unrestricted potential can be fairly arbitrary. Our approach is loosely connected with the probabilistic problem concerning the first-passage time of the Wiener process to a moving absorption boundary [23]- [25]. From a different perspective, our analysis has common points with the method of the so-called heat potentials as used in the theory of the parabolic partial differential equations [26,27]. We therefore give in the next section only the pivotal steps of the detailed mathematical treatment. Starting from the third section, we focus on the particular situation described in the first paragraph above.

Green function for the restricted diffusion
Consider the diffusion of a particle in the field of a space-homogenous and time-dependent force F(t), hence in the time-dependent potential V(x, t) = −xF(t). In the overdamped regime, the time-evolution of the probability density for the position of the diffusing particle is controlled by the Smoluchowski equation Here, v(t) is the time-dependent velocity of the drift motion, which is invoked by the external force, v(t) = F(t)/ , equals the particle mass times the viscous friction coefficient. The reciprocal value of the parameter is the particle mobility. The thermal-noise strength parameter D increases linearly with the absolute temperature, D = k B T/ . Let the initial condition be imposed at the time t , i.e. lim t→t p(x, t) = π(x). For the unrestricted diffusion, we assume the natural boundary conditions lim x→±∞ p(x, t) = 0, for any t t . The solution of the above Smoluchowski equation is well known. One possible treatment of the underlying time-ordered exponential makes use of the powerful algebraical methods, cf for example [19]. The result is This Green function yields the solution of equation (1) for the initial condition π(x) = δ(x − x ) imposed at time t . Qualitatively, it represents a gradually spreading Gauss curve whose centre moves in time. The momentary position of the maximum probability and at the same time, the mean particle position is controlled by the time protocol of the external drive. It is given by the The diffusive component of the probability current introduces an important function Let us consider this function as a function of the variable x. It exhibits an important property if x crosses the fixed coordinate x . Namely, it develops a singularity in the sense [27] S(x + , t; where δ(t − t ) is the Dirac distribution and is a positive infinitesimal quantity. Here and below, we tacitly assume the limit → 0 + . The limit will be performed in the expressions, where the presence of is unnecessary. The functions S(x ± , t; x , t ) always occur behind the time-integral sign. In this sense, we can treat them as operators acting on time-dependent functions. Owing to equation (4), these operators incorporate both an ordinary kernel and the Dirac distribution. More precisely, they act on an arbitrary time-dependent continuous function where S 0 (t, t ) is the ordinary kernel Obviously, in the limit t → t − , the analytical properties of the kernel S 0 (t, t ) are basically different from those of the kernels S(x ± , t; x , t ). Finally, for the sake of transparency of the following formulae, let us introduce the space derivative of the probability current, i.e. the function In closing the preparatory analysis, we wish to stress the following observation. The function S(x, t; x , t ) has been defined as the space-derivative of a solution of equation (1). Therefore, it represents per se also a solution of the Smoluchowski equation and will be as such used in the ansatz (8) and (9) for the general solution below. Let us now assume that an arbitrary fixed coordinate x = b divides the line into two half-lines. We assume the space-homogenous time-dependent forces F 1 (t) and F 2 (t) in the left half-line and the right half-line, respectively. We shall equip some of the above-defined functions with lower index to indicate, whether they will refer to the solution of the unrestricted diffusion problem with the force F 1 (t) (lower index '1') and with the force F 2 (t) (lower index '2'), respectively. Thus e.g. G 1 (x, t; x , t ) denotes the Green function for the unrestricted diffusion in the field of the space-homogenous force F 1 (t). Its explicit form follows from formula (2), if we use on the right-hand side v 1 (T ) = F 1 (t)/ instead of v(t). Let us further assume that the potential at the boundary exhibits a time-independent step of magnitude K > 0. The step is oriented downwards, if we move across the boundary from the left to the right. Finally, for simplicity, we assume that the initial condition π(x) = δ(x − x ) has been imposed at the time t = 0 and let x > b. Differently speaking, initially the particle is situated with probability one in the region to the right of the potential step.
We are interested in the Green function for the thus-defined diffusion problem. Let us designate it by U(x, t; x , 0), and let C(x, t; x , 0) be the corresponding probability current. If x b, the solution will be represented by the functions U 1 (x, t; x , 0), C 1 (x, t; x , 0) and similarly for the right half-line. We start with the ansatz where the functions N 1 (b, t , x ) and N 2 (b, t , x ) will be specified from the matching conditions. It is obvious that the functions (8) and (9) represent the solution of the Smoluchowski equation within the respective half-lines. Moreover, they satisfy the initial condition and the boundary conditions The assumed solution implies the probability current within the individual half-lines Let us now specify the matching conditions. The first matching condition concerns the probability current at the boundary. The current must be continuous. Hence, we require . This condition also guarantees the conservation of the total probability mass for the evolution governed by the Green function U(x, t; x , 0). This corollary emerges after the space-integration of the formulae (8), (9) and the time-derivative of the resulting expression. The second matching condition relates the values of the probability density at the points b ± . Having the potential step, the probability density suffers at x = b, a discontinuity. As a consequence of the first matching condition, one can prove Using the assumed solution (8)-(11), the matching conditions translate into a system of two integral equations for the unknown functions N 1 (b, t, x ) and N 2 (b, t, x ). Specifically, the integral equations read

6
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Generally speaking, the solution of this system yields through equations (8)-(11) the complete solution of the diffusion problem in question. However, in this paper, we are interested in the problem with a reflecting boundary. Within our formulation, the solution for the reflecting boundary can be achieved by simply invoking the limit of the infinitely high potential step. The limit K → +∞ implies ξ → 0 + and thereby also a significant simplification of the above integral equations. First of all, equation (12) reduces to a homogenous Volterra integral equation of the first kind. Its only solution is the trivial one. Hence, the function N 1 (b, t, x ) vanishes identically as function of the argument t. Thereupon, using equation (8), the probability density U 1 (x, t; x , 0) vanishes for any x b and t 0. Similarly, using equation (10), the probability current C 1 (x, t; x , 0) also vanishes, for any x b and t 0. The current across the boundary must be continuous and hence we immediately get C 2 (b + , t; x , 0) = 0. In other words, no trajectory of the emerging diffusion process can penetrate into the region to the left from the infinite step. Using the result that U 1 (x, t; x , 0) = 0 for any x < b, we get the conservation of the total probability mass within the right half-line. Summing up, the infinitely high potential step actually acts as the reflecting boundary.
Let us now turn to the second integral equation, equation (13). Taking into account that The equation can be further simplified. Using equation (7), the left-hand side reads In the last step, we have used G 2 (b + , t; b, t) = 0. We now time-integrate equation (14) and The equation will be solved in the following section.
In the rest of the present section, we wish to address the question of a possible meaning of the function N 2 (b, t, x ). In fact, we shall show that N 2 (b, t, x )/D represents the probability density at x = b for the problem with the reflecting boundary, i.e. we claim that The proof consists of two steps. First, we shall show that the Volterra integral equation of the first kind (16) is equivalent to the Volterra integral equation of the second kind, In order to prove the equivalence, we consider the above-formulated diffusion problem but without any potential step. Let H 1 (x, t; x , t) and H 2 (x, t; x , t) be the Green function for this 7 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT auxiliary problem. Once more, we assume the ansatz where M 1 (b, t, x ) and M 2 (b, t, x ) are two matching functions. Presently, both the probability density and the probability current must be continuous at x = b. We thus arrive at the system of two integral equations We now complete the definition of the auxiliary problem by a special assumption that the force F 1 (t) within the left half-line equals the force F 2 (t) to the right from the 'boundary'. Then, of course, the boundary is completely superfluous. Nevertheless, the above integral equations (20) and (21) must be still valid. Namely the ansatz for H 2 (x, t; x , 0), i.e. the right-hand side of equation (19), must give G 2 (x, t; x , 0). This is possible if and only if the function M 2 (b, t , x ) vanishes identically. Taking into account that R 1 (x, t; b, t ) is equal to R 2 (x, t; b, t ) and that this function is continuous across the 'boundary', the second matching equation (21) reads i.e. it is equivalent to equation (14). As for the first matching equation (20), we first apply On the left-hand side, we use the relation (4) and we arrive at the equation On the whole, analysing our auxiliary problem, we have proven the equivalence of the integral equations (22) and (24). At the same time, the equation (22) is equivalent to equation (14) and hence also to equation (16). Similarly, the Volterra integral equation of the second kind (24) is equivalent to the integral equation (17). Therefore the integral equations (16) and (17) are indeed equivalent.

Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
Let us return to the problem with the reflecting boundary. Making use of equation (9), the probability density at x = b + reads In the second step, we have used the integral equation (17). This finishes the proof of the statement concerning the meaning of the matching function N 2 (b, t, x ).
In closing this section, we wish to introduce some notational simplifications and review the main formulae to be used in the following. The quantities with the lower index '1' which have referred to the left half-line will no more occur in the paper. Therefore, we can suppress the lower index '2' for the quantities, which have referred to the right half-line. For example, we shall use U(x, t; x , 0) instead of U 2 (x, t; x , 0) and G(x, t; x , t ) instead of G 2 (x, t; x , t ). Further, starting from this point, the boundary will be always placed at the origin, i.e. we take b = 0. Hence, there could be no confusion, if we use the designation N(t, x ) instead of N 2 (b = 0, t, x ). Using this agreement, the construction of the Green function U(x, t; x , 0) for the problem with time-dependent force and with the reflecting boundary involves three basic steps. Firstly, given the detailed time protocol of the force F(t), we calculate the Green function G(x, t; x , t ) for the spatially unrestricted diffusion. The result is given in equation (2). Secondly, we have to solve the Volterra integral equation of the first kind The solution of this equation describes through N(t, x )/D = U(0, t; x , 0), the probability density at the origin for the problem with the reflecting boundary. At the same time, making use of the relation the functions N(t, x ) and G(x, t; 0, t ) yield the complete solution of the problem. The probability current C(x, t; x , 0) vanishes at x = 0. The function U(x, t; x , 0) is properly normalized, i.e. we have ∞ 0 U(x, t; x , 0) dx = 1 for any t and for x > b.

Oscillating force with time-independent component
We are primarily interested in the response to the oscillating external force. A physically relevant situation occurs, if the symmetrically oscillating component F 1 sin(ωt) is superimposed on the time-independent negative force F 0 < 0, which pushes the particle against the reflecting boundary at the origin. Altogether, the force to be discussed in the rest of the paper will be taken as The three parameters occurring in this formula together with the diffusion constant D form the full description of our setting.

9
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Probability density at the boundary
Before embarking on further discussion, let us collect some intuitively expected features of the solution of equation (16). If F 1 = 0 or if ω = 0, we have the standard diffusion problem with a time-independent drift and a reflecting boundary at the origin. In this case, the kernel in the integral equation (16) depends only on the time-difference and the Laplace-transformation method [30] readily yields the explicit result Here, v 0 = F 0 / is the drift velocity and erfc(·) is the (complementary) error integral [31]. The last formula is valid for any F 0 . If F 0 0, the time-asymptotic value of the solution is zero because the particle gradually escapes far to the right from the boundary. On the other hand, for any F 0 < 0, one has lim t→∞ N(t, x )/D = |v 0 |/D. In this case, the force F 0 pushes the particle against the reflecting boundary and an arbitrary initial condition relaxes to the equilibrium state π eq (x) Returning to the general case, assume again the negative static force F 0 < 0, but let us have now F 1 = 0, ω = 0. Then the integral equation (16) represents a fairly non-trivial problem. There are two sources of difficulties. Firstly, the integrated function displays a weak singularity at the upper integration limit. Secondly and more importantly, the kernel in (16) is not of the convolution type. This observation rules out the application of the Laplace transformation method. In order to get some further insight, let us exceptionally introduce the reduced variables φ 0 = v 0 /(2Da), φ 1 = v 1 /(2Da), τ = a 2 Dt, = ω/(a 2 D), ξ = ax /2, where 1/a is the unit length. After several purely algebraical steps the final form of the integral equation (16) reads After solving the equation for X(ξ , τ), the probability density at the origin will follow from the transformation U(0, t; x , 0) = N(t, x )/D = aX(ξ , τ). Equations of the above form have been already investigated in the mathematical literature. As for their numerical solution, we have implemented the so-called Product Integration method [28,29], which includes a special treatment of the weak singularity in equation (30). Figure 1 illustrates the solution for several combinations of the parameters involved. In all four panels, the full curve describes the density U(0, t; x , 0) as obtained by the numerical solution of the above integral equation. The dashed curve is the density at the boundary for the corresponding problem without modulation, i.e. the function (29)  half-period of every such cycle the force F 1 sin(ωt) has pushed the particle to the right from the boundary. This implies a delayed decrease of the probability density at the boundary. If there is a phase shift in the input force, the response would also acquire a corresponding phase shift but its amplitude would remain unchanged. Comparison of the panels (a) and (b) (or (c) and (d), respectively) demonstrates the influence of increasing the amplitude of the external force. The oscillations are more pronounced and the response exhibits strong nonlinear features. Generally speaking, a highly non-equilibrium stationary regime occurs whenever the ratio F 1 /|F 0 | is greater than one, the temperature is low and the frequency of the force modulation is small.

Mean position
Let us now focus on the mean position of the particle. If F 0 < 0 and F 1 = 0, the probability density relaxes to the above-mentioned exponential form π eq (x). Thereupon, the particle mean position will eventually approach the value µ eq = D/|v 0 |. Turning to the case with the modulated force, we start with the general solution as given in equation (27). The mean position is defined as the first moment of the probability density. Hence, we have The first term on the right-hand side, denoted as µ 1 (t, x ), is simply the (conditional) mean position for the unrestricted diffusion, under the condition of x 0. This term exhibits damped nonlinear oscillations and finally vanishes. Actually, considering the unrestricted diffusion with negative constant force component as in equation (28), the particle will eventually glide towards −∞. This reasoning can be supported by the explicit formula for the function µ 1 (t, x ). But it is not very interesting and we omit it. Let us designate the second term on the right-hand side of equation (31) as µ 2 (t, x ) and the curly-bracketed expression in it as M(t, t ). Using equation (3) and carrying out the space integration, we get Thus the second term in equation (31) is the t -time integral of the product M(t, t )N(t , x ). Notice however, that the integral again does not have a convolution structure. Further, it still incorporates (through the function N(t , x )) transitory effects. We are of course primarily interested in the stationary regime. However, it does not seem possible to derive an appropriate asymptotic expansion of the function µ 2 (t, x ), which would allow an explicit elimination of the transitory effects. In any case, the stationary component of the non-convolution integral in question emerges in the long time limit. Formally, we can write As indicated, this function does not depend on the initial condition x . It exhibits oscillations which however are not centred around the equilibrium mean position µ eq = D/|v 0 |. Introducing the time-averaged mean position in the stationary regime as the dynamical shift of the mean position equals the difference µ a − µ eq . The dynamical shift is always positive. It always increases with the amplitude F 1 of the modulated force and it always decreases with the frequency of the force and with the temperature. Figure 2 shows the time-dependent mean position for the same set of parameters as used in figure 1.

Internal energy
Up to now, we have discussed kinetic quantities. They are directly related to the time-and spacedependence of the probability density. Recently, it was revealed that the kinetic behaviour can qualitatively differ from the analysis of the underlying energetics [15].
Let us again start with the equilibrium situation. As already mentioned, in the timeindependent potential V(x) = −xF 0 = x|F 0 |, the system approaches the equilibrium state described by the probability density π eq (x) = (|v 0 |/D) exp[−x|v 0 |/D]. The equilibrium internal energy is simply the spatial integral of the product of the functions V(x) and π eq (x). The result is E eq = D = k B T . Hence, the internal energy does not depend on the slope of the potential and it linearly increases with the temperature. Similarly, in equilibrium, the (configurational) Helmholtz free energy is F eq = k B T log[|F 0 |/(k B T )] and the entropy follows from the relation F eq = E eq − TS eq . In our specific setting, the changes of the Helmholtz free energy are purely of entropic origin, since the equilibrium internal energy E eq is independent on the 'width'of the potential.
Returning to the time-dependent force, the argumentation of the stochastic energetics [14] is based on the Langevin equation. The position of the particle is represented by the random variable X(t). If the particle occurs at an arbitrary fixed instant at the position X(t), it has the potential energy V [X(t), t], where V(x, t) = −xF(t) is the time-dependent potential of the problem. From this perspective, having solved the Smoluchowski equation (1), we know the probability density for the random variable X(t)-it is just the Green function (27) from the preceding section. The mean internal energy is then the probabilistic mean value of the random variable V [X(t), t], i.e. we have where the mean position µ(t, x ) is given in the formula (31). Generally speaking, the internal energy exhibits all transitory effects included in the time evolution of the mean position. In the stationary regime, it is again a periodic function of time with the fundamental period ω/(2π). The oscillations express the combined effect of both the periodically modulated heat flow to the bath and the periodic exchange of the work done on the system by the external agent. We are interested in the time-averaged internal energy in the stationary regime, i.e. in the time-independent quantity This quantity is always greater than the equilibrium mean energy E eq = D . Figure 3 shows the difference E a − E eq as the function of the diffusion constant D and the amplitude of the external force F 1 .

Heat and work
The above applied methods of the stochastic energetics incorporate also the separate calculation of both the heat and the work. Generally speaking, the heat (≡ the dissipated energy) can be identified as the 'work' done by the particle on the heat bath. It arises if and only if the particle moves, i.e. it is inevitable connected with the probability density current. More precisely, within our setting, the heat released to the heat bath during the time interval [0, t] is given as where I(t, x ) = ∞ 0 dx J(x, t; x , 0) is the integrated current. On the other hand, the external agent does work on the system by lifting the potential V(x, t), while the position of the particle is virtually fixed. Thus the work done at a given instant depends on the momentary position of the particle. Averaging over the possible positions, the work done be the external field during the time interval [0, t] reads where µ(t, x ) is again the mean position of the particle (31). Similarly to the preceding subsection, we now define the average work done by the external agent during one fundamental This time-independent quantity must be always positive. Otherwise, the system in contact with the single heat bath at one fixed temperature would produce positive work during the cyclic process in question. Numerically, the average work equals the area inside the hysteresis curve which represents the parametric plot of the oscillating force F 1 (t) versus the mean coordinate µ a (t). Figure 4 presents the time-averaged work W a as the function of the temperature and the amplitude F 1 . The parameters used in figure 4 are the same as those in figure 3.

Conclusions
One of the aims of this paper was to explore the flow of energy from its external source through the observed system to its final destination, which is the heat bath.
In the stationary regime, the system exhibits periodic changes of its state. Therefore its internal energy, as every state function, is also periodic. This of course says nothing about the period-averaged value of the internal energy. However, we have shown that the time-averaged energy content in the stationary non-equilibrium regime, when the heat is being released at a fixed bath-temperature, is always higher than the equilibrium internal energy corresponding to the same temperature. In other words, the system can act as an energy-flow channel only by having a permanently increased time-averaged internal energy.
Having periodic changes of the internal energy, the work done on the system during one period must be equal to the heat dissipated during the period. However, their behaviour during an infinitesimal time interval within the period is quite different. The heat released to the heat bath is always positive. Actually, at any given instant, the force F(t) and the motion of the particle have the same direction. Hence, the function which is integrated in the last expression in equation (37) is always non-negative.
Contrary to this, the work done on the system during an infinitesimal time interval can be both positive and negative. In order to be specific, consider again the force F 0 + F 1 sin(ωt) with the negative component F 0 (the amplitude F 1 is always positive). Let W (i) a , i = 1, . . . , 4, denote the work done by the external field during the ith quarter-period of the force modulation. During the first quarter-period the slope of the potential decreases and the particle does a work, irrespective to its momentary position. Thus we have W (1) a < 0. Nevertheless, the farther the particle from the boundary, the bigger the work done by it during the fixed time interval. Within the second and the third quarter-period, the slope of the potential increases and the work is done by the external agent. Hence, we have W (2) a > 0 and W (3) a > 0. However W (2) a is bigger than