Phase space reduction of the one-dimensional Fokker-Planck (Kramers) equation

Phase space reduction of the one-dimensional Fokker-Planck (Kramers) equation

Abstract

A pointlike particle of finite mass , moving in a one-dimensional viscous environment and biased by a spatially dependent force, is considered. We present a rigorous mapping of the Fokker-Planck equation, which determines evolution of the particle density in phase space, onto the spatial coordinate . The result is the Smoluchowski equation, valid in the overdamped limit, , with a series of corrections expanded in powers of . They are determined unambiguously within the recurrence mapping procedure. The method and the results are interpreted on the simplest model with no field and on the damped harmonic oscillator.

pacs:
05.40.Jc, 87.10.Ed

I I. Introduction

The Brownian motion of a particle in a confined system represents an essential model used in description of stochastic transport through quasi one-dimensional (1D) systems, e.g. channels in nanomaterials, pores or fibers in biological structures. In a 1D system, the trajectory of a Brownian particle is described by the Langevin equation

(1.1)

Here, denotes mass of the particle, represents the driving potential, is an effective friction coefficient and is the stochastic force, satisfying the usual conditions on averaged values, and ; is the temperature and the Boltzmann constant. The corresponding phase space density of the particle satisfies the Fokker-Planck (FP), or Kramers (kinetic) equation,

(1.3)

where is the spatial coordinate, denotes its velocity and is the inverse temperature.

Solutions of the Langevin equation, as well as the corresponding kinetic equation, are studied over almost a century (1); (2); (3). Still, motion of a particle in confined geometries represents usually a complicated problem, requiring next reductions of the used description. Due to simplicity, and often for relevance, mainly the overdamped limit is studied. The mass dependent term in Eq. (1.1), , is then considered negligible and the particle’s spatial density is governed by the Smoluchowski equation,

(1.4)

containing no information about the mass of the particle; denotes the diffusion constant. Instead of the full phase space, one works with only the spatial coordinate . Of course, the question is the price for this simplification, as well as possibility of also including properly inertia of the massive particles in the reduced (real space only) description of the Brownian dynamics. Recent studies showed its importance in understanding rectification of the transport in ratchets (4), or its influence on the effective diffusion coefficient (5); (6); (7) in a narrow channel. For our purpose to demonstrate the phase space reduction, we will deal here with only the 1D FP (Kramers) equation (1.3).

The eq. (1.3) is reducible to the Smoluchowski eq. (1.4) in the limit . The reduction procedure (8); (9); (10) is based on an instant thermalization of the particle’s velocity after any move in the direction in the case of an infinitely small mass. The situation resembles derivation of the Fick-Jacobs (FJ) equation (11); (12),

(1.5)

appearing as the result of the dimensional reduction of the diffusion equation in a 2D channel with varying cross section , onto the longitudinal coordinate , if equilibration in the transverse () direction is instant. The function denotes the linear (1D) density of the particle. Of course, as in reduction of diffusion to the FJ equation, that of (1.3) to (1.4) as is a singular limit, and must be handled with care, but from the viewpoint of the reduction of (1.1), no such caveat is needed.

Recently, an exact mapping procedure was proposed (13); (14), enabling us to also derive the corrections to the FJ equation (1.5), which are necessary, if the transverse equilibration is not instant. The key was to introduce anisotropy of the diffusion constant in the diffusion equation,

(1.6)

governing the 2D (spatial) density . For , the infinitely fast transverse diffusion immediately flattens the profile of . Then integration of Eq. (1.6) over the cross section, together with the reflecting boundary conditions satisfied at the hard walls, results in the FJ eq. (1.5). In the case of a slower transverse diffusion, the mapping procedure generates a series of corrections to the FJ equation controlled by , developing the profile of , which is already curved depending on the flux and geometry of the channel.

The procedure was extended to mapping of diffusion in an external field , e.g. for diffusion in a channel with soft walls (15), where the particle is kept near the axis by the parabolic potential, ; represents the varying stiffness of the walls. The equation to be mapped onto the coordinate is the 2D Smoluchowski equation,

(1.7)
(1.8)

Integrating over and applying the mapping scheme gives the mapped 1D equation in an extended Smoluchowski form,

(1.9)

governing the 1D density , where stands for an effective potential and is the correction operator sought as an expansion in the small parameter . The potential and the coefficients of depend on and both are fixed unambiguously within the recurrence mapping procedure.

The central idea of this paper is a conjecture that the FP eq. (1.3) can be reduced to an extended Smoluchowski-like 1D form, governing the spatial density , after integration over and applying a similar mapping scheme. The velocity thus represents a ”transverse” coordinate instead of , with the mass playing the role of the small parameter . Then the Smoluchowski equation (1.4) should be obtained in the limit (8); (9). Performing the recurrence procedure, a series of corrections to this equation in powers of would be derived. Then the final mapped equation would also respect inertia of the Brownian particles, although working only in real space.

Use of the mapping procedure developed for diffusion (13); (14); (15); (16) is not straightforward; the left-hand side operator of the FP equation (1.3) has a different structure than that of the diffusion (1.6) or Smoluchowski equation (1.7). Still, there is a way to apply the general scheme of the mapping in this case and to perform reduction of the phase space onto the real space in the way described above. Presentation of this algorithm is the primary aim of this study.

The result of the mapping of the FP equation (1.3) is again an equation of the form (1.9) with and replaced by . In the limit of stationary flow, this equation can be simplified by subsequent reduction of the operator to a function , a spatially dependent effective diffusion coefficient (12); (17),

(1.10)

The leading term of is proportional to and the whole series of corrections to the Smoluchowski equation can be summed up, giving

(1.11)

with , if the higher derivatives of are neglected.

In the following section, we analyze how to apply the mapping scheme for reduction of the phase space in the FP equation (1.3). In Sect. III, our considerations are verified on an exactly solvable model, the FP equation with no field, . Analysis of this example helps us to construct the recurrence scheme for calculation of the series of corrections to the zeroth order approximation, the Smoluchowski eq. (1.4), in the small parameter . Finally, the complete mapping procedure for an arbitrary potential is presented in Sect. IV. The mapped equation of type (1.9), as well as the formula (1.11) is derived, and checked on the damped linear harmonic oscillator.

Ii II. Preliminary considerations

The key points of the mapping procedure, as formulated for diffusion (13); (14), are recalled in this Section. Based on physical considerations, we adjust the general scheme of the mapping for dimensional reduction of the FP equation (1.3) to this situation.

The mapping procedure represents a consistent transition from a fine grain to a coarse grain description of some evolution process. The process is described in details by some partial differential equation (PDE), governing the density of particles in the full space, defined by the coordinates . The dimensional reduction projects this equation onto another PDE, which governs the density in the reduced space of the coordinate . The coordinate is one of the coordinates of the full space, , and the mapping accomplishes integration over the transverse coordinates . In the case of the FP equation (1.3), phase space represents the full space and the dimensional reduction integrates over the ”transverse” coordinate, the velocity . Hence we have the defining relation between the densities and :

(2.1)

The phase-space density is expected to be near the thermodynamic equilibrium and so approximately proportional to , which provides convergence of the integral in Eq. (2.1).

Then the first step of the mapping is also integration of Eq. (1.3) over . If completed, we get

(2.2)

the other terms are zero due to in the limit . This equation represents nothing but mass conservation in the reduced space,

(2.3)

where the 1D flux is defined by the relation

(2.4)

In contrast to diffusion, where is fixed, here the flux is a function formally independent of the density , so we also need the evolution equation for this quantity. After integration of Eq. (1.3) multiplied by , we obtain

(2.5)
(2.6)

[after some algebra and applying the definitions (2.1) and (2.4). This step recalls the Grad’s method of moments (18); (19). For the 1D Kramers equation (1.3), taking only a couple of the zeroth () and the first () order moment of the phase space density is satisfactory for generating the selfconsistent system of the mapped (real space) equations (2.3) and (2.5).]

The next key point of the mapping algorithm is that of expressing the full-space density using the 1D density and also the flux in this case. This relation enables us to complete integrations in Eq. (2.5) and get the evolution equation for , together with the 1D mass conservation, in closed form. The initial task is to find the zeroth order approximation, valid in the limit . Our first proposal for such a relation between and the reduced space quantities and is based on the following physical construction:

For an infinitely small mass of the particle, the stochastic force thermalizes its velocity almost immediately after any move along the spatial coordinate . Similar to the transverse equilibration of the 2D density of a particle diffusing in a narrow channel with biasing transverse force (15); (20); (21); (22), one could try the formula with separated Boltzmann factor in the fast relaxing ”direction” , . It is easy to check that it does not work here; the flux becomes zero according to Eq. (2.4). To prevent this failure, let us suppose that the distribution in is shifted by the local mean (macroscopic) velocity , depending on the local flux, . Then we have

(2.7)

Retaining only these two terms in the square brackets and replacing by the flux , one gets

(2.8)

which will be taken as the sought zeroth order relation between and the reduced space quantities, and .

This heuristic formula will be verified later by the exact mapping algorithm. Still, one can check immediately that the relation (2.8) represents correctly a kind of backward mapping of the 1D (spatial) functions and onto the phase space densities ; if substituted for in the defining relations (2.1) and (2.4), we obtain identities. Applying Eq. (2.8) to Eq. (2.5), the integrals over can be completed and the result,

(2.9)

together with the mass conservation, Eq. (2.3), forms a closed couple of PDE, governing the mapped quantities and .

In the limit , the first term in Eq. (2.9), , becomes negligible and the equation expresses the zeroth order relation between the flux and the density ,

(2.10)

If combined with the mass conservation, Eq. (2.3), we get the Smoluchowski equation (1.4); represents the diffusion constant.

The calculation presented shows that the Smoluchowski equation (1.4) is related to the FP equation (1.3) in the same way as the Fick-Jacobs equation (11) to the diffusion equation valid in a narrow 2D channel. Both mapped equations describe an asymptotic behavior of the full-space density infinitely rapidly equilibrating in the transverse direction; the velocity plays the role of the transverse coordinate for the FP equation. Our considerations indicate that the mass of the particle, , becomes the small parameter, controlling the series of corrections to the Smoluchowski equation in the case when the transverse equilibration is not infinitely fast.

Following the scheme of the mapping procedure (13); (14), the next point is that of searching for the true relation between the phase space density and the 1D quantities and , replacing the heuristic formula (2.8), valid for nonzero . Without losing generality, it can be written in the form

(2.11)
(2.12)

If the operators and (with implicit) are expandable in , one can substitute for in the FP equation (1.3) and fix the coefficients of these operators to satisfy this equation in each order of , similar to the mapping of diffusion. Then, using the relation of backward mapping (2.11) in Eq. (2.5) gives the expansion of the evolution equation for and finally, in combination with mass conservation (2.3), the sought series of corrections to the Smoluchowski equation (1.4) in terms of the finite mass .

To verify whether this scheme is viable, we analyze the exactly solvable case with in the next Section.

Iii III. Exactly solvable model

The exact solution of the FP equation (1.3) with no potential, , is presented in this Section. We demonstrate the mapping on the example of the phase space density evolving from the initial density with thermalized velocity . The mapped equation, as well as the form of the operators and in Eq. (2.11) can be found explicitly in this case. This analysis will direct us in construction of the recurrence mapping scheme in Sect. IV.

The case is exactly solvable (23), the Green’s function of the FP equation (1.3),

(3.1)
(3.2)

can be derived explicitly (see Appendix A),

(3.5)

if expressed in the scaled coordinates,

(3.6)
(3.7)
(3.8)

and . If the thermalized particle is inserted at time with a spatial distribution ,

(3.9)

evolution of the density is given by the formula

(3.11)
(3.13)

, the abbreviations

(3.14)
(3.15)

are used and the integration over runs over the whole (unspecified) 1D spatial domain.

Then the spatial (1D) density and the flux are integrated directly according to Eqs. (2.1) and (2.4),

(3.16)
(3.18)

It is easy to check that the mass conservation (2.3), in the scaled coordinates, is satisfied. The quantity plays the role of a ”stretched” time (24). For short times, , and . The formulas (3.11) and (3.16) describe correctly behavior of the Newtonian particles in this limit. The mapping procedure, as outlined in the previous Section, requires us to study asymptotic behavior in the opposite limit, .

For large times, , the stretched time becomes and the formulas (3.16) represent the general solution of the diffusion equation, as expected according to the previous Section. Now it is necessary to verify that the mass can serve as a small parameter controlling the series of corrections to the diffusion equation and its solution.

It may look problematic at first glance, because contains , representing essential singularity of the variable . Then the formulas (3.16) [and similarly Eq. (3.11)] are not expandable in . Nevertheless, this property is still consistent with the general scheme of the mapping, as analyzed in Ref. (14).

The dimensional reduction, as demonstrated on anisotropic diffusion in a narrow channel (14), also reduces the Hilbert space of the full-space problem. Let us denote the spatial operator of the evolution equation, i.e. for anisotropic 2D diffusion; the eigenvalues and the eigenfunctions are given by the equation

(3.19)

supplemented by proper boundary conditions at the walls of the channel. The parameter of anisotropy splits the spectrum into two parts, the low-lying states, whose eigenvalues remain finite for and the transients with diverging . Then the exact 2D density evolves as

(3.20)

the constants are given by the initial condition and the summation runs over the whole spectrum. The transients contribute to the sum by the terms proportional to , where are finite in the limit . The result is a formula containing the essential singularity in the parameter near zero, similar to Eq. (3.11) with singular for .

The mapping procedure reduces the full Hilbert space of all onto the space defined only by the low-lying states . If the 1D density is integrated from Eq. (3.20) and mapped backward onto the full Hilbert space (by some operator ), the transients will be canceled; the sum (3.20) after the mapping there and back runs only over the low-lying states . The terms retained involve no essential singularity in ; the formula for considered in the mapping procedure represents the regular part of the exact 2D density with respect to the parameter near zero. This reduction is natural for the zero-th order (Fick-Jacobs) approximation, as the transients decay infinitely fast due to their infinite eigenvalues . Nevertheless, the mapping based on fixing the series of corrections expanded in can work only with the regular part of the 2D density.

Correspondingly, the formulas (3.11) and (3.16) are exact, including the contributions of the transients, which are represented by the singular terms . The mapping requires us to analyze only the regular parts,

(3.22)
(3.24)

and

(3.26)

written in the unscaled coordinates, obtained after taking only the regular parts of and (3.14), and , in Eqs. (3.16) and (3.11). Of course, the formulas are applied for , when the transients vanish. We omit writing the subscript ”reg” in the following text.

In comparison with the overdamped limit , evolution of the spatial (1D) density and the flux is only corrected by a time shift, in the formulas (3.22). The Gaussian distribution is retarded by the time , corresponding to the mean time necessary for losing information about the original velocity. The value of the shift is constant in and , so the density (3.22) still obeys the diffusion equation,

(3.27)

and . There is no correction of the zero-th order mapped equation due to the finite mass of the particle. Let us stress that the case of is extremely simple, similar to the mapping of the 2D diffusion in a flat narrow channel, also giving no corrections to the Fick-Jacobs approximation.

Nevertheless, the relation of the backward mapping, generating the phase-space density from the mapped quantities and , is not quite trivial. In the formula for , Eq. (3.26), the time is shifted by and the displacement by with respect to the distribution of a massless particle. The shortest way to construct the relation of backward mapping is by applying the shift operators in and on , compensating the different shifts of time and displacement in Eqs. (3.22) and (3.26),

(3.28)

Due to Eq. (3.27), the operator in the exponent can be replaced by . After expanding the shift operator and using , we arrive at the relation of the form (2.11),

(3.30)

the explicit formulas for and to be substituted in Eq. (2.11) in the case are

(3.31)
(3.32)

Both operators are expandable in , (); their zero-th order coefficients equal unity, consistent with the heuristic formula (2.8). Also applying the relation (3.30) in the definitions (2.1) and (2.4) gives identity.

The final step is that of verifying the evolution equation for , Eq. (2.5). Two integrals are to be completed with expressed by the backward mapping, Eq. (3.30),

(3.33)
(3.34)

which are now valid exactly. If substituted in Eq. (2.5), and the equation (3.27) is applied, we get

(3.35)

after simple algebra. This equation validates the relation for nonzero as well and thus, if combined with the mass conservation (2.3), also the diffusion equation (3.27) without any corrections. If compared with the calculation of Eq. (2.10), the 1-st order term , neglected in the previous Section, is compensated here by other 1-st order term coming from the exact relation (3.30), appearing in the integrals (3.33).

The FP equation (1.3) with zero potential is too simple to give nonzero corrections to the diffusion equation (3.27). Nevertheless, it helped us to understand the structure of the mapping. It shows that the scheme suggested in the previous section is viable. The relation of the backward mapping has the form of Eq. (2.11), at least for , and the operators and can be expanded in , or ,

(3.36)

Integration in Eqs. (3.33) indicates that the coefficients and should be sought dependent up on the scaled velocity rather than (compare to Ref. (25)); otherwise each term in Eq. (3.36) would contribute in several succeeding orders in the integrals (3.33). The mixing of orders would hinder us in constructing the recurrence scheme generating corrections to Eq. (2.5). In the notation of Eq. (3.36), and

(3.37)
(3.38)
(3.39)
(3.40)

valid for according to Eqs. (3.31), will be used for testing the results of the recurrence procedure in the next Section.

Iv IV. Mapping procedure

We now finish the construction of the mapping procedure, outlined in the Section II, for an arbitrary analytic potential . Supposing the backward mapping of the form (2.11) with the operators and expanded in () according to Eqs. (3.36), we find recurrence relations fixing the coefficients and . Completing the integrals in Eq. (2.5), we obtain a series of corrections to the zero-th order relation . Combined with mass conservation, Eq. (2.3), it gives the Smoluchowski equation corrected due to nonzero mass of the particle.

The essential relation determining the operators and is the FP equation (1.3), which has to be satisfied for any solution of the reduced problem, the density and the flux , after their backward mapping (2.11) onto the full-dimensional Hilbert space. If the expansion in of both operators, Eq. (3.36), is supposed, we have

(4.1)
(4.2)

after introducing the scaled velocity in Eq. (1.3). The factor is replaced by in the following calculations. Thus half-integer powers of appear in Eq. (4.1) (25). As this equation has to be satisfied for any , we can split it for clarity into two relations: the first one, including only the integer powers of ,

(4.3)
(4.4)

and the second one, collecting the half-integer powers,

(4.5)
(4.6)

Notice that Eqs. (4.3) and (4.5) do not violate parity of and in . If used for construction of the recurrence relations between the coefficients, all they have to have the same parity as ; hence and . This symmetry enables us to find the normalization (or identity) conditions for and . The backward mapped , Eq. (2.11), with the operators , expanded in , Eq. (3.36), substituted in the definitions (2.1) and (2.4) has to give identities for any , and . Thus we obtain

(4.7)
(4.8)

The operators and are supposed not to depend on time, so the time derivative commutes with them and acts directly on or in Eqs. (4.3), (4.5). To derive the operators unambiguously, using only spatial derivatives, we express from the mass conservation, Eq. (2.3). However, the time derivative of cannot be expressed in a similar way from Eq. (2.5), because is not the leading term there. If the backward mapping, Eqs. (2.11) and (3.36), is applied to the integrals of Eq. (2.5), we get

where the operators are given by

(4.9)

the right-hand side integral of Eq. (2.5) results in

(4.10)

after integrating by parts and using the normalization relation (4.8). Then, instead of the evolution equation for , Eq. (2.5) has to be understood as an expression relating and ,

(4.12)

and the flux , as well as its time derivative, is expressed using according to this equation. The term acts now on as the 1-st order correction in after completing inversion and commutation with the spatial operators in the square brackets of Eq. (4.12). Then is systematically replaced by , contributing to the next corrections in the higher orders of . The result is a formula for expressed by some purely spatial operator acting on ,

(4.13)

where the operators are related to , Eq. (4.9). Using this form, we are able to write the final mapped equation explicitly after combination with Eq. (2.3),

(4.14)
(4.15)

which is the Smoluchowski equation (1.4) in the zero-th order, extended by a series of mass dependent corrections in . The operators are introduced to simplify notation in the following calculations.

Now the operators can be expressed explicitly using . Expanding in and applying Eq. (4.14), we have

which is to be used in the operator equation,

(4.16)
(4.17)

obtained after comparison of Eqs. (4.12) and (4.13). Expanding in and comparing the coefficients of the same powers of , we get a sequence of relations fixing ,