How to solve Fokker-Planck equation treating mixed eigenvalue spectrum?

An analogy of the Fokker-Planck equation (FPE) with the Schr\"odinger equation allows us to use quantum mechanics technique to find the analytical solution of the FPE in a number of cases. However, previous studies have been limited to the Schr\"odinger potential with a discrete eigenvalue spectrum. Here, we will show how this approach can be also applied to a mixed eigenvalue spectrum with bounded and free states. We solve the FPE with boundaries located at x=\pm L/2 and take the limit L\rightarrow\infty, considering the examples with constant Schr\"{o}dinger potential and with P\"{o}schl-Teller potential. An oversimplified approach was proposed earlier by M.T. Araujo and E. Drigo Filho. A detailed investigation of the two examples shows that the correct solution, obtained in this paper, is consistent with the expected Fokker-Planck dynamics.


Introduction
The one-dimensional Fokker-Planck equation (FPE) for the probability density p(x, t), depending on variable x and time t, assumes the generic form [1][2][3][4][5][6][7] ∂p(x, t) ∂t = − ∂ ∂x {f (x, t)p(x, t)} + ∂ 2 ∂x 2 D(x, t) 2 p(x, t) . The FPE provides a very useful tool for modeling a wide variety of stochastic phenomena arising in physics, chemistry, biology, finance, traffic flow, etc. [3][4][5][6]. Given the importance of the Fokker-Planck equation, different analytical and numerical methods have been proposed for its solution. As it is well known, the stationary solution of FPE can be given in a closed form if the condition of detailed balance holds. The study of the time-dependent solution is a much more complicated problem. The FPE (1.1) with a general time-dependence and a special x-dependence of the drift and diffusion coefficients has been studied analytically in [7] using Lie algebra. This method is applicable when the Fokker-Planck equation has certain algebraic structure, allowing to apply the Lie algebra and the Wei-Norman theorem. Generally, there are few exactly solvable cases. A simple such example is a system with constant diffusion coefficient and harmonic interaction of the form f (x) = −dV (x)/dx with harmonic potential V (x) ∼ x 2 . The case with double-well potential is already quite non-trivial and requires a numerical approach [8].
The known relation between the Fokker-Planck equation and the Schrödinger equation can also be used. This approach allows us to apply the well known methods of quantum mechanics. In particular, analytical solutions can be found in the cases, where the eigenvalues and eigenfunctions for the considered Schrödinger potential are known. For a general Schrödinger potential, numerical treatments used in quantum mechanics, such as the Crank-Nicolson time propagation with implicit Numerov's method for second order derivatives [9], are very useful. To apply it to Schrödinger-type equation, we just need to replace the real time step ∆t by an imaginary time step ∆t → −i∆t. In quantum mechanics this is called imaginary time propagation and is used for calculation of ground states and also excited states. The analytical studies of mapping the FPE to Schrödinger equation has been up to now restricted to a treatment of discrete eigenstates. An attempt has been made in [10] to extend this approach to potentials with mixed (discrete and continuous) eigenvalue spectrum. However, we have found a basic error in this treatment, indicated explicitly at the end Sec. 4.3.
The aim of our work is to show how the problem with mixed eigenvalue spectrum can be treated correctly. We will show this in two examples: one with constant Schrödinger potential and another -with Pöschl-Teller potential. The same example has been incorrectly treated in [10]. To avoid any confusion one has to note that the Pöschl-Teller potential is called Rosen-Morse potential in [10].

Solution of FPE with constant diffusion coefficient
We start our consideration with the one-dimensional Fokker-Planck equation (1.1) for the probability density distribution p(x, t), depending on variable x and time t. Here f (x) is the nonlinear force and D is the diffusion coefficient, which is now assumed to be constant. We consider natural boundary conditions and take the most frequently used initial condition in the form of the δ-function. This FPE (2.1) can be transformed into an equation of Schrödinger type (see Sec. 2.2). Unfortunately, the well known relation (see Eq. (2.25)), derived for the discrete eigenvalue spectrum, cannot be applied if this equation has continuous or mixed eigenvalue spectrum. To overcome this problem, we follow a properly corrected treatment of [10]. Namely, we solve the FPE with boundaries located at x = ±L/2 and then take the limit L → ∞ (see Sec. 2.3). This approach is applied in quantum mechanics to describe unbounded states. To keep closer touch with quantum mechanics, here we will use the boundary conditions p(x = ±L/2, t) = 0, further called absorbing boundaries.

The stationary solution
The stationary solution p st (x) is the long-time limit of p(x, t) at t → ∞, which follows from the equation The force f (x) can be expressed in terms of the potential V (x) via f (x) = −dV (x)/dx. It yields (2.5) ?????-2 Due to the natural boundary conditions, we have zero flux Thus, we have which yields the stationary solution has meaning of the unnormalized stationary solution only in case of natural boundaries and N is the normalization constant This function Y (x) is further used to construct the time-dependent solution.

The time-dependent solution with discrete eigenvalues
Here we derive the time-dependent solution, starting with the transformation p(x, t) → q(x, t) defined by This transformation removes the first derivative in the original Fokker-Planck equation and generates the equation of Schrödinger type for the function q(x, t), i. e., is the so-called Schrödinger potential. In the case of discrete eigenvalues, we apply the superposition ansatz After inserting (2.15) into (2.13), we get the eigenvalue problem for eigenfunctions ψ n (x) and eigenvalues λ n ≥ 0 with time-dependent coefficients a n (t) given by a n (t) = a n (0) exp (−λ n t) . (2.17) ?????-3 According to this, Eq. (2.15) can be written as and satisfy the closure condition (completeness relation) Eq. (2.16) can be written as a Schrödinger-type eigenvalue equation with Hermitian Hamilton operator H: The coefficients a n (0) in (2.18) are calculated using the initial condition According to (2.18), this relation can be written as In the following, we multiply both sides of this equation by ψ n (x) and integrate over x from −∞ to +∞. Taking into account (2.19), it yields the up to now unknown coefficients The final result of this calculation reads Note that this method can also be used for other boundary conditions. The solution in the general form of (2.25) is well known from older studies, e. g., [11] and can be found in many textbooks, e. g., [3,4].

The time-dependent solution with mixed eigenvalue spectrum
Consider now the problem with two absorbing boundaries located at x = ±L/2 instead of the natural boundary conditions. In this case we have a discrete eigenvalue spectrum, and Eq. (2.25) can be used (with summation over only those eigenfunctions which satisfy boundary conditions in a box of length L) to calculate the probability distribution p L (x, t), i. e., where λ n,L are eigenvalues and ψ n,L (x) are corresponding eigenfunctions, which fulfill the boundary conditions. Let us split this infinite sum into two parts: for λ n,L < λ con and λ n,L ≥ λ con , where λ con ?????-4 is the smallest continuum eigenvalue in the case with natural boundaries. This eigenvalue spectrum is shown schematically in Fig. 1, where the value of λ con is indicated by a horizontal dotted line, the eigenvalues λ n,L < λ con -by solid lines and the eigenvalues λ n,L ≥ λ con -by dashed lines.
Let M (L) be maximal value of n for which λ n,L < λ con and k n−M (L),L = 2 D (λ n,L − λ con ) for n > M (L) and ψ con k n−M (L),L (x) = ψ n,L (x) for n > M (L). Hence, we have The solution with natural boundaries is the limit case L → ∞ where the normalization constant N is given by The infinite sum can be split in two parts: one with odd m and the other -with even m. If the Schrödinger potential is symmetric, then one of these two parts contains only odd eigenfunctions ψ o k (x), whereas the other part -only even eigenfunctionsψ e k (x). In the limit L → ∞, these two sums can be represented by corresponding integrals, yielding where g e (k) = lim This representation is useful if the eigenvalues and eigenfunctions are known.

The analytical solution of FPE with constant force
Let us consider a constant force term. In this case the Fokker-Planck equation (2.1) reads This is a drift-diffusion problem for the potential V (x) = −v drift x normalized to V (x = 0) = 0.
No stationary solution exists for this problem, because the normalization constant N in Eq. (2.11) diverges in this case. Nevertheless, the transformation (2.12) p( can be used here to obtain an equation of Schrödinger type (2.13) with constant Schrödinger potential The corresponding to (2.21) stationary Schrödinger-type equation reads where n = 0, 1, 2, . . . and k n,L = π L (n + 1). (3.7) According to (3.6)-(3.7), we have from (2.33) and (2.34) Taking into account that holds, we obtain from Eq. (2.32) the expression [cos(kx) cos(kx 0 ) + sin(kx) sin(kx 0 )] . we obtain after a simplification the well known result which describes a moving and broadening Gaussian profile.

Fokker-Planck dynamics with Pöschl-Teller potential
Here we consider as a particular example the force with some positive constants b and α. This corresponds to the diffusion problem in the potential normalized to V (x = 0) = 0. Fig. 2 shows that this potential is actually a smoothed version of the V-shaped potential. The corresponding Schrödinger potential in this case is . If we compare it (see Eq. (4.4) and Fig. 3) with the well known Pöschl-Teller potential then we see that Eq. (4.3) represents the shifted by b 2 2D Pöschl-Teller potential with V 0 = b 2 2D + bα 2 . As we can see from Fig. 3, the Pöschl-Teller potential gives mixed (discrete and continuous) eigenvalue spectrum, therefore Eq. (2.25) cannot be directly applied to solve the FPE. We have to use (2.32).

Bounded solutions for Pöschl-Teller potential
The Eq. The eigenvalues can be calculated from the following equation [12] λ n = −(l − n) 2 , for n < N ; n ∈ N . (4.7) Note that at least one bounded state withλ 0 = −l 2 always exists forl > 0, which corresponds to λ 0 = 0. The bounded eigenfunctions are known [12] where F denotes hypergeometric function, which can be represented by Gaussian hypergeometric series The normalization constants are
As we see, the eigenfunctions are rather complicated in general case. The expressions become essentially simpler for integer values ofl. Therefore, without loosing the general idea, we will show the solutions of the Fokker-Planck equation forl = 1 andl = 2. The unbounded eigenfunctions (4.13) and (4.14) arē Inserting these relations and also λ con =l 2 α 2 D/2 (following fromλ con = 2λcon Dα 2 −l 2 = 0) into (2.32), we finally obtain the time-dependent solution of the Fokker-Planck equation

The solution of FPE for
If the initial condition is given by x 0 = 0, then ψ õ k (0) = 0 and ψ ẽ k (0) = 1 hold, which allows us to obtain a simpler expression (4.31) The solution for parameters b = 2, D = 2 and α = 1, corresponding tol = 1, with the initial location of the delta-peak at x 0 = 5 is shown in Fig. 4 for different time moments t. As we  The stationary solution is practically reached at t = 10. This behavior is expected from the driftdiffusion dynamics. For small times t → 0, we have a delta-peak located at x = x 0 in accordance with the given initial condition (2.3). For comparison, the "general solution" of [10] does not satisfy this initial condition, as a result of a wrong construction, where the contribution of bounded states is simply summed up with a Gaussian probability density profile (calculated with an error). The latter corresponds to unbounded states for zero Schrödinger potential at L → ∞, as it is evident from (3.13) and (3.3) at v drift = 0. Therefore also the result appears to be correct only at t → ∞ when the Gaussian part vanishes. It is clear that the whole set of eigenfunctions must be calculated self-consistently for the given potential to obtain a correct and meaningful result, since only in this case the completeness relation (2.20) holds and all different eigenfunctions are orthogonal.
????? -11 Thus, the basic error of [10] is that some of the eigenfunctions are calculated for zero Schrödinger potential in [10], whereas all of them must be calculated for the true Schrödinger potential.

The solution of FPE for Pöschl-Teller potential with parameterl = 2
Forl = 2 (which implies b = 2αD) we have two bounded states with eigenvaluesλ 0 = −4 and λ 1 = −1. The corresponding eigenfunctions are (4.32) The unbounded eigenfunctions arē By adding again two absorbing boundaries atx = ±L/2, we have discrete values ofk, i. e., kL ,m for even functions andκL ,m for odd functions. In the limitL → ∞, we obtain again the classical infinite-square-well relations for eigenstates: The normalization constants in this case are  The solution for parameters b = 4, D = 2 and α = 1, corresponding tol = 2, with the initial condition given by x 0 = 5 is shown in Fig. 5 for different time moments t. The evolution of the probability distribution is very similar to that one shown in Fig. 4 forl = 1, with the only essential difference that the dynamics is faster and the distribution is somewhat narrower because of a deeper potential well.

Conclusions
Using the analogy of the Fokker-Planck equation with the Schrödinger equation, it has been shown how the time-dependent solution can be constructed in the case of mixed eigenvalue spectrum with free and bounded states. The method is based on the idea of introducing two absorbing boundaries at x = ±L/2, considering the limit L → ∞ afterwards. Although this idea is similar to the one proposed earlier in [10], it is obvious that the problem is quite non-trivial, so that the oversimplified (erroneous) approach of [10] cannot be used -see discussion at the end of Sec. 4.3. Analytical solutions have been found and analyzed in two examples when the Schrödinger potential is constant (constant force) and a shifted Pöschl-Teller potential. For the latter potential, the analytical solutions have been compared with the results of the Crank-Nicolson numerical integration method, and the agreement within an error of 10 −7 has been found. The time evolution of the calculated probability distribution in these examples is consistent with the usual drift-diffusion dynamics.