Solution of the Fokker-Planck equation with a logarithmic potential and mixed eigenvalue spectrum

TL;DR

This paper solves a Fokker-Planck equation with a logarithmic potential related to Bessel processes, revealing a mixed eigenvalue spectrum and introducing a new normalization technique for continuous eigenfunctions, with applications in physics, biology, and finance.

Contribution

It presents a novel method to evaluate the normalization of continuous eigenfunctions in Fokker-Planck equations with mixed spectra, applied to Bessel-like processes with negative drift.

Findings

Eigenfunction normalization technique based on asymptotic behavior

Solution of Fokker-Planck with mixed eigenvalue spectrum

Comparison with other analytical and numerical methods

Abstract

Motivated by a problem in climate dynamics, we investigate the solution of a Bessel-like process with negative constant drift, described by a Fokker-Planck equation with a potential V(x) = - [b \ln(x) + a\, x], for b>0 and a<0. The problem belongs to a family of Fokker-Planck equations with logarithmic potentials closely related to the Bessel process, that has been extensively studied for its applications in physics, biology and finance. The Bessel-like process we consider can be solved by seeking solutions through an expansion into a complete set of eigenfunctions. The associated imaginary-time Schroedinger equation exhibits a mix of discrete and continuous eigenvalue spectra, corresponding to the quantum Coulomb potential describing the bound states of the hydrogen atom. We present a technique to evaluate the normalization factor of the continuous spectrum of eigenfunctions that relies solely upon their asymptotic behavior. We demonstrate the technique by solving the Brownian motion problem and the Bessel process both with a negative constant drift. We conclude with a comparison with other analytical methods and with numerical solutions.

Figures

Equations

$${\partial }_t\rho (x,t)=-{\partial }_x\left[f(x)\rho (x,t)\right]+\frac{1}{2}{\sigma }^2{\partial }_x^2\rho (x,t),$$

$$\rho (x,t)=\frac{1}{(2\pi {\sigma }^{2}t{)}^{\frac{1}{2}}}{\text{e}}^{-\frac{x^2}{2{\sigma }^{2}t}},$$

$${\rho }_{st}(x)=\frac{\mathcal{N}}{{\sigma }^2}{\text{e}}^{-\frac{2}{{\sigma }^2}V(x)},$$

$${\partial }_t\rho (x,t)=-{\partial }_x\left[-\frac{b}{x}\rho (x,t)\right]+\frac{1}{2}{\sigma }^2{\partial }_x^2\rho (x,t),$$

$$\rho (x,t)=x^{\left(\frac{1}{2}-\frac{b}{{\sigma }^2}\right)}{x}_0^{\frac{1}{2}+\frac{b}{{\sigma }^2}}\frac{1}{{\sigma }^{2}t}{\text{e}}^{-\frac{x^2+x_0^2}{2{\sigma }^{2}t}}{I}_{\frac{1}{2}+\frac{b}{{\sigma }^2}}\left(\frac{x{x}_0}{{\sigma }^{2}t}\right),$$

$$V(x)=\frac{b}{2}\ln (1+x^2),\text{and hence}f(x)=-\frac{bx}{1+x^2},$$

$$V(x)=b\ln (x)+\frac{a}{2}{x}^2,\text{and hence}f(x)=-\frac{b}{x}-ax,$$

$${\partial }_t\rho (x,t)=-{\partial }_x\left[\left(\frac{b}{x}+a\right)\rho (x,t)\right]+\frac{1}{2}{\sigma }^2{\partial }_x^2\rho (x,t),$$

$${\int }_{x_1}^{x_2}\rho (x,t)dx=\frac{A}{R}$$

$${\partial }_t\rho (x,t)=\mathcal{L}\rho (x,t);\mathcal{L}=-{\partial }_{x}f(x,t)+\frac{1}{2}{\partial }_x^{2}g(x,t{)}^2,$$

$${\partial }_t\rho (x,t)=\mathcal{L}\rho (x,t);\mathcal{L}=-{\partial }_{x}f(x)+\frac{1}{2}{\sigma }^2{\partial }_x^2,$$

$$dx=f(x)dt+\sigma dW(t),$$

$$\rho (x,t)={\int }_{\Omega }d{x}_0\rho (x,t|x_0,t_0){\rho }_{in}(x_0,t_0),$$

$$\rho (x,t|x_0,t_0)={\text{e}}^{\mathcal{L}(t-t_0)}\delta (x-x_0).$$

$${\rho }_{st}(x)=\frac{\mathcal{N}}{{\sigma }^2}{\text{e}}^{-\frac{2}{{\sigma }^2}V(x)},$$

$${\int }_{\Omega }dx{\rho }_{st}(x)=1,$$

$$\rho (x,t|x_0,t_0)=\sum _{\lambda \in \Lambda }{\text{e}}^{-\lambda (t-t_0)}{\varphi }_{\lambda }(x_0){\varphi }_{\lambda }(x),$$

$$\mathcal{L}{\varphi }_{\lambda }(x)=-\lambda {\varphi }_{\lambda }(x),$$

$$\rho (x,t)={\int }_{\Omega }d{x}_0\rho (x,t|x_0,t_0){\rho }_{in}(x_0,t_0)=\sum _{\lambda \in \Lambda }{\alpha }_{\lambda }(t_0){\text{e}}^{-\lambda (t-t_0)}{\varphi }_{\lambda }(x),$$

$${\alpha }_{\lambda }(t_0)={\int }_{\Omega }d{x}_0{\rho }_{in}(x_0,t_0){\varphi }_{\lambda }(x_0).$$

$$\mathcal{L}=\frac{{\sigma }^2}{2}{\partial }_x{\text{e}}^{-\frac{2}{{\sigma }^2}V(x)}{\partial }_x{\text{e}}^{\frac{2}{{\sigma }^2}V(x)}$$

$${\int }_{\Omega }dx{\varphi }_1(x)\left\{\mathcal{L}{\varphi }_2(x)\right\}\ne {\int }_{\Omega }dx{\varphi }_2(x)\left\{\mathcal{L}{\varphi }_1(x)\right\}.$$

$$\tilde{\mathcal{L}}={\text{e}}^{\frac{2}{{\sigma }^2}V(x)}\mathcal{L},$$

$${\varphi }_{\lambda }(x)={\text{e}}^{-\frac{V(x)}{{\sigma }^2}}{\psi }_{\lambda }(x)\propto \sqrt{{\rho }_{st}(x)}{\psi }_{\lambda }(x).$$

$$\text{scalebox}0.8\underset{t\to t_0}{\lim }\rho (x,t|x_0,t_0)=\delta (x-x_0)=\sum _{\lambda \in \Lambda }{\psi }_{\lambda }(x_0){\psi }_{\lambda }(x)={\text{e}}^{\frac{V(x)}{{\sigma }^2}+\frac{V(x_0)}{{\sigma }^2}}\sum _{\lambda \in \Lambda }{\varphi }_{\lambda }(x_0){\varphi }_{\lambda }(x),$$

$$\delta (x-x_0)={\text{e}}^{-\frac{V(x)}{{\sigma }^2}+\frac{V(x_0)}{{\sigma }^2}}\delta (x-x_0)={\text{e}}^{+\frac{V(x)}{{\sigma }^2}-\frac{V(x_0)}{{\sigma }^2}}\delta (x-x_0)$$

$$={\text{e}}^{\frac{2}{{\sigma }^2}V(x)}\sum _{\lambda \in \Lambda }{\varphi }_{\lambda }(x_0){\varphi }_{\lambda }(x)={\text{e}}^{\frac{2}{{\sigma }^2}V(x_0)}\sum _{\lambda \in \Lambda }{\varphi }_{\lambda }(x_0){\varphi }_{\lambda }(x).$$

$$\rho (x,t)={\int }_{\Omega }d{x}_0\sum _{\lambda \in \Lambda }{\text{e}}^{-\lambda (t-t_0)}\sqrt{\frac{{\rho }_{st}(x)}{{\rho }_{st}(x_0)}}{\psi }_{\lambda }(x_0){\psi }_{\lambda }(x){\rho }_{in}(x_0,t_0).$$

$$-\lambda {\psi }_{\lambda }(x)=\tilde{\mathcal{L}}{\psi }_{\lambda }(x)=\left[\frac{{\sigma }^2}{2}{\partial }_x^2-{\mathcal{V}}_s(x)\right]{\psi }_{\lambda }(x),$$

$${\mathcal{V}}_s(x)\equiv -\frac{1}{2}\left[d_x^{2}V(x)-\frac{1}{{\sigma }^2}(d_{x}V(x){)}^2\right]=\frac{1}{2}\left[d_{x}f(x)+\frac{1}{{\sigma }^2}f(x{)}^2\right].$$

$${\int }_{\Omega }dx{\psi }_n(x){\psi }_m(x)=0,m\ne n,$$

Full text

Fokker-Planck equation, Bessel process, Schrödinger equation, Logarithmic potential, mixed eigenvalue spectrum, $\delta$-function normalization.

††preprint: preprint

I Introduction

Many complex systems can be investigated by capturing the effect of a large number of degrees of freedom by adding a noise term to a deterministic evolution equation. This may be done either phenomenologically or by implementing a stochastic reduction procedure on a more complex set of equations (e.g. franzke2005, ; MW2017, ). Statistical information is governed by the probability density function $\rho(x,t)$, which satisfies the Fokker-Planck equation

[TABLE]

a second order parabolic partial differential equation wherein $f(x)$ describes the deterministic dynamics and $\sigma$ is the constant noise strength. (Note, for parsimony of notation, we write the operator $\partial^{2}/\partial{x^{2}}$ as $\partial_{x}^{2}$.) For example, when a one-dimensional solely diffusive process governs $\rho(x,t)$, there is no deterministic drift force ($f(x)=0$), and the Fokker-Planck equation is a diffusion equation with a Gaussian solution;

[TABLE]

for initial data $\rho(x,t=0)=\delta(x)$, with natural boundary conditions $\rho(x=0,t)=\rho(x=\infty,t)=0$111Here natural, absorbing and reflecting boundaries refer to boundaries where the probability density vanishes (sufficiently fast to insure normalization), vanishes with finite flux, or has zero flux, respectively. Typically, natural boundary conditions are imposed at $\pm\infty$. See also Feller FELL1952 . The long-time behavior of Eq. (1) for $f(x)\neq 0$ is given by the Boltzmann distribution in the case of thermal systems, defined as

[TABLE]

where $\mathcal{N}$ is a normalization constant and we can consider $V(x)=-\int^{x}f(y)\,dy$ to be a potential in which, say, a Brownian particle evolves and hence the connotation of the deterministic drift force; $f(x)=-d_{x}V(x)$.

A large class of systems are described by a logarithmic potential (sometimes referred to as an entropy term), viz., $V(x)=b\,\ln(x)$, with $b$ constant. Depending on the sign of $b$, the central force associated with this potential is either attractive or repulsive, and systems governed by this potential are generally referred to as Bessel processes, characterized by the Fokker-Planck equation

[TABLE]

which has the following exact solution

[TABLE]

in which $I_{\frac{1}{2}+\frac{b}{\sigma^{2}}}$ is the modified Bessel function of the first kind, and $x_{0}=x(t=0)$. Due to the singularity at the origin, Eq.(4) is defined on the half-line [0,$\infty$), the origin satisfies a reflecting boundary condition and there is a natural boundary condition $\rho(x=\infty,t)=0$ at infinity.

The Bessel process exhibits several interesting properties. For example, as described in (refs. BRWL2000, ; MBGB2010, ), by integrating over angle variables and letting $b\to b-(d-1)\,\sigma^{2}/2$, the equation for the radial component in $d$-dimensions is equivalent to the one-dimensional form given by Eq. (4). Thus, a particle in a logarithmic potential in a space dimension $d$ will have the same radial distribution function as a free-particle in $\tilde{d}=d-2\,b/\sigma^{2}$ dimensions, where $\tilde{d}$ can be negative and non-integer. Pitman and Yor PIYO1981 ; PIYO1982 studied the relationships between Brownian motion and Bessel processes, and squared Bessel processes are connected to local times (or occupation times) in Brownian motion through a space-time transformations and the application of Ray-Knight theorems (e.g., YORE1991, ; YOMA1999, ). Hence, squared Bessel processes play a major role in the study of maps between Brownian motion in different dimensions WILL1970 ; PITM1975 .

The importance of Bessel diffusion processes is generally two-fold. Firstly, it belongs to a restricted class of exactly solvable models, along with Brownian motion and the Ornstein-Uhlenbeck process. Secondly, it is widely applicable in physics, chemistry, biology and finance. For example, in physics, it describes (a) the effective long-range interaction between two probe particles in a one-dimensional driven fluid, such as solute particles in a solvent LMSL2005 ; (b) the vortex-antivortex annihilation in the 2D XY Heisenberg model subject to thermal fluctuations at temperatures below the Kosterlitz-Thouless transition $T_{KT}$ BRWL2000 ; (c) an alternative technique to quantize non-Abelian gauge field theories, avoiding the drawbacks of the Fadeev-Popov method HHSA1983 222Note that in Ref. HHSA1983 the Bessel process is mistakenly referred to as the Rayleigh process.. In finance, the property that the Bessel process is equivalent to exponential Brownian motion with linear drift in time (or geometric Brownian motion) makes it a suitable tool for the evaluation of Asian options, which are contracts based on the average underlying price and are thus considered less expensive than standard options, due to the lower variance of the averages GEYO1993 ; GJYO2003 .

As summarized nicely by Dechant et al., (DLBK2011, ), because the Bessel process with a purely logarithmic potential $V(x)=b\,\ln|x|$ is singular at the origin, then the steady state given by the Boltzmann distribution of Eq. (3) is non-integrable. This is because $\rho_{st}(x)\propto\text{e}^{-\frac{2}{\sigma^{2}}\,V(x)}=x^{-2\,b/\sigma^{2}}$ clearly diverges at large $x$ (small $x$) depending on whether $2\,b/\sigma^{2}<1$ ($2\,b/\sigma^{2}>1$). However, it is more common that the potential of real systems only behaves like a pure logarithm asymptotically. Therefore, it is more appropriate to talk about “Bessel-like processes”. For example, in the settings of (a) and (b) discussed above, there is an effective characteristic scale $a$ such that at large $x$, the potential behaves as $V(x)\approx b\,\ln(x/a)$. In other settings, such as in the case of the diffusive spreading of momenta in two-level atoms in optical lattices DGRT2006 , the potential may also exhibit a more complex logarithmic dependence, such as

[TABLE]

where $x$ is the momentum.

Because of the slow power law spatial decay, few moments of the stationary distribution will diverge. Therefore, complete information regarding the stationary distribution is not provided by the time-independent solution $\rho_{st}$, but instead requires taking the infinite-time limit of the time-dependent moments DLBK2011 . The latter has been proved by Lutz to be related to a breaking in the ergodicity for systems with power-law distributions LUTP2004 . By adding regular terms to the potential, the existence of the stationary solution can always be guaranteed. One example is the Rayleigh process RAYL1880 ; RAYL1905 ; SNTR1963 ; GNRS1986 , given by

[TABLE]

which is considered to be a generalized Bessel process, or a Bessel process with a linear drift or radial Ornstein-Uhlenbeck process, as it generalizes the radial projection of a $d$-dimensional Ornstein-Uhlenbeck equation ORUH1930 . The Rayleigh process is widely used in finance because the parabolic term corresponds to the Cox-Ingersoll-Ross (CIR) process CIRO1985 . The CIR process describes the evolution of financial variables that are strictly positive and empirically known to exhibit mean reversion, in which there is a convergence to the average value in the long time limit GJYO2003 .

Here we focus on another generalization of the Bessel process wherein the drift is characterized by the potential $V(x)=-b\ln(x)-a\,x$, and is thus associated with the Fokker-Planck equation,

[TABLE]

which arises in applications in biology, queuing theory, finance and climate dynamics. In biology, it models the dynamics of DNA bubbles forming when the hydrogen bonds of the base pairs in the Watson-Crick double-helix structure break due to thermal fluctuations FMDD2007 . The continuous zipping and unzipping of the double-strand, referred to as “DNA breathing”, allows chemicals and proteins to bind to reactive sites of the bases that they would not otherwise be able to access FMDD2007 . Thus, when $x$ in Eq.(8) refers to the dynamics of the bubble length, this is referred to as the Poland-Scheraga model. In queueing theory, Cofmann et al. (CPRE1998, ) contrast the relevance of the Bessel process and reflective Brownian motion–a random walk with a reflective singularity at the origin. The former (latter) case describes the distribution of the unfinished work in the heavy-traffic limit for a two queue system served by a single server that switches between them instantaneously (over a finite time). In finance, this process has been applied to the evaluation of Asian options VALI2004 . Linetsky LINE2004 uses a non-linear transformation and a measure change of Eq. (8) to find a generalization of the CIR process that provides a better fit of empirical data. In climate dynamics, Toppaladoddi and Wettlaufer have derived a Fokker-Planck-like equation for the distribution of Arctic sea ice thickness that, in geophysically relevant limits, can be described by a Bessel process with constant drift, and have found both stationary and time-dependent solutions to match satellite data STJW2015 ; STJW2017 . The concept of the sea ice thickness distribution was introduced by Thorndike et al. Thorndike:1975 , and is defined by considering a region with area $R$ that is sufficiently large to contain a range ice of different thicknesses. Then the integral

[TABLE]

gives the fraction of that area ($A/R$) that contains ice of thicknesses between $x_{1}$ and $x_{2}$. Then, for a simple treatment of the thermal growth process, the spatio-temporal evolution of $\rho(x,t)$, subject to wind, thermal and mechanical forcing, is governed by Eq. (8). A fully coupled approach to the climatological evolution of the ice thickness distribution uses a thermodynamic only nonlinear nonautonomous deterministic backbone EWNT2009 and thus requires a numerical approach. However, a stochastic perturbation theory has been developed that provides numerical tests of the Langevin equations emerging from a range of generalizations of the geophysical models MWJMP2013 ; STJW2015 ; MWJC2017 .

Whilst the “bare” Bessel process has been solved exactly, the Bessel process with constant drift poses significant challenges in this regard. Progress can be made by solving the Fokker-Planck equation in terms of either Fourier or Legendre transforms and numerically transforming back LNLE2014 , or by expanding the solution directly into a complete set of orthonormal eigenfunctions and solving the associated Sturm-Liouville eigenvalue problem. For an infinite domain, the presence of a constant drift term leads to an eigenvalue spectrum that is a mix of a discrete and continuous set of eigenvalues, the latter of which are associated with eigenfunctions that oscillate at infinity, posing obvious complications. Indeed, oscillating solutions must be normalized to a $\delta$-function, and thus to unity in the sense of distributions. The process described by Eq. (8) has been solved with several techniques. Linetsky gave a spectral representation both by solving the equation in real space in a finite box and then sending the box size to infinity VALI2004 , and by solving the equation in complex space in the full domain $[0,\infty)$ LINE2004 (see TITC1962 for an introduction to the technique). We note that the associated Sturm-Liouville equation strongly resembles the Schrödinger equation for the Coulomb interaction in quantum mechanics, describing the discrete eigenstates of the hydrogen atom. Here, motivated by the climate dynamics problem described above, we employ a different technique to solve this problem. The method evaluates the normalization of oscillatory eigenfunctions on the infinite half-line without the use of complex variables, and solely requires knowledge of their asymptotic oscillatory behavior (heuristically like the method of stationary phase, but clearly without appeal to the Riemann-Lebesgue Lemma). The paper is organized as follows. In §II we introduce the spectral expansion of the Fokker-Planck equation, and describe how to normalize oscillating solutions. In §III, we use this method to solve the problem of Brownian motion with constant negative drift. Finally, in §IV we solve the Bessel process with negative constant negative drift, and compare the results with other analytic techniques and numerical solutions.

II Spectral Expansion of the Fokker-Planck equation

The general form of the Fokker-Planck equation in the Itô formulation is

[TABLE]

where $\mathcal{L}$ is the Fokker-Planck operator, and $f(x,t)$ and $g(x,t)$ are the drift and diffusion functions respectively. For given initial data, Eq. (10) satisfies conditions at the boundary of the domain $\Omega=[x_{in},x_{out}]$. We will consider the simpler case of additive noise, where $g(x)\equiv\sigma^{2}$ defines the constant noise strength, and time-independent constant drift $f(x,t)\equiv f(x)$, and hence the Fokker-Planck equation becomes

[TABLE]

corresponding to which is the autonomous Langevin equation

[TABLE]

where $dW(t)$ is a Wiener process. The probability density can be written as

[TABLE]

where $\rho_{in}(x_{0},t_{0})$ is an initial condition at time $t_{0}$ and position $x_{0}$. The transition density satisfies Eq. (11) with initial condition \raisebox{2.15277pt}{\scalebox{0.8}{\displaystyle\lim_{t\to t_{0}};}}\rho(x,t|x_{0},t_{0})=\delta(x-x_{0}), and thus formally

[TABLE]

Thus, as described in §I, the stationary solution of Eq. (11) is

[TABLE]

where $\mathcal{N}$ is a normalization constant, viz.,

[TABLE]

and we can consider $V(x)=-\int^{x}f(y)\,dy$ to be a potential associated with the deterministic drift force; $f(x)=-d_{x}V(x)$. The existence of a renormalizable stationary solution is guaranteed by requiring that the potential diverges at both boundaries, and hence for example a particle is strictly confined within the domain.

II.1 Eigenfunction expansion

The Fokker-Planck equation can be solved exactly just in few cases. Generally, a solution can be obtained in real space by expanding the transition density into a set of eigenfunctions and solving the associated eigenvalue problem RFPE1934 . We expand the transition density as

[TABLE]

where the $\phi_{\lambda}$ are the set of eigenfunctions of $\mathcal{L}$ satisfying the eigenvalue problem

[TABLE]

and $\Lambda$ is the eigenvalue spectrum. Inserting the expansion (17) into the density (13) we find

[TABLE]

where the coefficients $\alpha_{\lambda}(t_{0})$ are the projections of the initial condition on the $\lambda^{\text{th}}$ function, viz.,

[TABLE]

The calculation of the spectrum can be simplified by rewriting Eq. (18) in terms of a self-adjoint operator. The Sturm-Liouville theorem ensures that such operators have a complete basis of orthonormal eigenfunctions and a real spectrum. However, the operator $\mathcal{L}$ in Eq. (18) is not self-adjoint, which can be seen by writing $\mathcal{L}$ in terms of the potential $V(x)$ as

[TABLE]

and checking that for two test functions $\phi_{1}$ and $\phi_{2}$ the following inequality holds

[TABLE]

Nevertheless, we can construct a self-adjoint operator associated with the Fokker-Planck operator as

[TABLE]

which satisfies Eq. (18) and possesses a complete orthonormal set of eigenfunctions $\{\psi_{\lambda}\}$ with non-negative real eigenvalues $\lambda\geq 0$. The two operators have the same spectrum, and their eigenfunctions are, up to normalization, related by

[TABLE]

The completeness of the basis insures that the initial condition of the transition density satisfies

[TABLE]

which, by adding unity, in the sense of distributions, can be rewritten as

[TABLE]

Substituting (26) into (14), and then into (19) using (24), we find

[TABLE]

II.2 The Schrödinger-like equation

With the aid of Eq. (24) we can reparametrize the Fokker-Planck equation (11) into the following Sturm-Liouville problem

[TABLE]

which resembles an imaginary-time Schrödinger equation with potential

[TABLE]

Again, the Sturm-Liouville theorem guarantees that the set of eigenfunctions, $\psi_{\lambda}(x)$, is orthonormal and complete, and all eigenvalues are real. In the case of a finite domain, the spectrum is an infinite set of discrete eigenvalues. However, if we take the limit of an infinite domain size the eigenvalue spectrum can either (a) remain a discrete set $\Lambda_{d}$, (b) become a continuum $\Lambda_{c}$, or (c) be a mix of the two $\Lambda=\Lambda_{d}+\Lambda_{c}$. Intuitively, as shown in Fig. 1, the type of spectrum in the infinite domain case can be understood by depicting eigenvalues as lines of constant Schrödinger potential, Eq. (30), or Schrödinger isopleths. Thus, a region of the spectrum wherein isopleths are bound at both extremes of the potential exhibits a discrete set of eigenvalues. Hence, since the operator $\mathcal{\tilde{L}}$ is self-adjoint, the spectrum in this region can be solved simply by imposing the boundary conditions on the generic solutions $\psi_{\lambda}(x)=\text{C}\,u_{\lambda}(x)$, where C is an integration constant. The discrete set of solutions satisfying the boundary conditions will occur at values $\lambda_{n}$, $n=0,1,2,\cdots$, with $\lambda_{n=0}=0$ corresponding to the stationary solution, when it exists. Each discrete solution oscillates in the domain with $n$ zeros and is orthogonal to any solution with $m\neq n$, namely,

[TABLE]

where $n\equiv\lambda_{n}$ and $m\equiv\lambda_{m}$ are two eigenvalues. The integration constant is determined by normalization of the solution as

[TABLE]

to satisfy the orthonormality relation

[TABLE]

where $\delta_{nm}$ is the Kronecker delta. In case A of Fig.1 all isopleths are bound, the entire spectrum is discrete, and the solution of the Fokker-Planck equation is given by Eq. (28).

When solved in a finite domain, a region of the spectrum wherein the potential bounds one or none of a line’s extremes also exhibits a discrete set of eigenvalues. The spectrum is thus solved using the approach just described. However, as shown in case C of Fig.1, when taking the boundaries to infinity, the set of lines becomes dense and tends towards a continuum. Therefore, except special cases, there is a continuous set of solutions, $\psi^{con}$, that oscillate asympotically at large $x$ with a finite constant amplitude. Oscillating functions cannot be normalized in the sense of Eq. (33), as the integral does not exist. That said, so long as they are real, they can be normalized in the sense of distributions to a Dirac $\delta$-function as

[TABLE]

where $k_{1}$ and $k_{2}$ are two eigenvalues in the continuous set. When the entire spectrum is a continuum, the sum in Eq. (28) passes to an integral and the transition density is

[TABLE]

Finally, at one (or both) of the boundaries approaching infinity, the potential may asymptotically approach a critical value $\lambda_{cr}$ that demarcates the continuous from the discrete region of the spectrum, as shown in panel B of Fig.1. Hence, we expect either a finite or an infinite set of discrete eigenfunctions $\psi_{\lambda_{n}}^{dis}$, with eigenvalues $\lambda_{n}$, $n=0,1,2,\cdots$ for $\lambda<\lambda_{cr}$, and a continuous set, $\psi_{k}^{con}$, with eigenvalues $k\in[0,\infty)$, where $\lambda=\lambda_{cr}+k$. Thus, solution for the transition density in this case is

[TABLE]

where $N$ is the highest discrete eigenvalue for a finite discrete set.

II.3 Normalization of oscillating functions

As mentioned in §II.2, the oscillatory solutions cannot be normalized in the sense of Eq. (32). Nevertheless, the normalization provided by Eq. (34) does not provide an intuitive way to calculate the normalization constant. One approach is to solve the Sturm-Liouville problem in a finite domain of size $L$. The continuum of oscillating functions is then reduced to a discrete set of functions satisfying the boundary conditions of the problem. Then, the functions can be normalized in the sense of Eq. (32), after which we can take the limit $L\to\infty$ as in BKMH2013 ; BLIS1975 . Here we adapt the approach introduced in the 1920s by Fues FUES1926 ; BSQM1957 to normalize quantum mechanical oscillatory eigenfunctions, which provides further insight into the structure of the Sturm-Liouville solutions. We extend the expansion of the solution shown in Eq. (19) to the case of a mixed spectrum as

[TABLE]

in which $N$ and $\lambda_{cr}$ are the same as in Eq.(II.2), $a_{k}(t_{0})$ is the projection of the initial condition on the $k$-th eigenfunction, and $\phi(x)\propto\sqrt{\rho_{st}(x)}\,\psi(x)$. We then split the integral into a series of infinitesimal integrals viz.,

[TABLE]

where the $\Delta_{n}k$, $n=1,2,\cdots,\infty$, form an infinite discrete set of eigenvalue shells. Here, $k_{n}$ denotes the eigenvalue of the $n^{\text{th}}$ shell such that \raisebox{2.15277pt}{\scalebox{0.8}{\displaystyle\lim_{\Delta_{n}k\to 0};}}k_{n}=k, and the functions $\Psi_{n}(x)$ are the eigendifferentials, defined as the integrals of eigenfunctions over the $n^{\text{th}}$ shell;

[TABLE]

Although a single eigenfunction $\psi_{k}^{con}(x)$ oscillates asymptotically in its argument, a linear combination of such functions does not, and hence the eigendifferentials decay at large $x$. Therefore, eigendifferentials play the same role in the continuum that the discrete eigenfunctions play in their spectral region. Consequently, discrete eigenfunctions are orthogonal to eigendifferentials,

[TABLE]

for all $n$ and $m$. Thus, we can generalize the orthonormality definition to

[TABLE]

Determination of the factor $\Delta_{n}k$ in the orthonormalization rule for $\Psi_{n}(x)$ is given in Appendix A. The projection $a_{k}(t_{0})$ in Eq. (37) can then be defined as

[TABLE]

Now, by using the eigendifferential definition of Eq. (39) in the third of Eqs. (41), we obtain for $n=m$

[TABLE]

where $\{k_{1},k_{2}\}\in\Delta_{n}k$. The prefactor $(\Delta_{n}k)^{-1}$ insures that the above integral holds for any arbitrary shell size so long as

[TABLE]

where $u(k,x)$ is the general solution of the Sturm-Liouville equation and $\text{C}(k)$ the normalization. Therefore, Eq. (44) generalizes Eq. (32) to oscillatory solutions.

Next we show that the normalization constant in Eq. (44) depends solely on the amplitude of the asymptotic oscillations of $\psi_{k}^{con}(x)$. Let us take the Schrödinger-like equation (29) for two different eigenvalues $k_{1}$ and $k_{2}$, multiply them by $\psi^{con}_{k_{2}}$ and $\psi^{con}_{k_{1}}$ respectively, and then subtract them, which gives

[TABLE]

Now, using Eq. (45) in Eq. (44), and dividing by $(k_{1}-k_{2})$, leads to

[TABLE]

We integrate by parts in a finite domain with the boundary terms defined at infinity, and consider two domains; (i) \raisebox{2.15277pt}{\scalebox{0.8}{\displaystyle\lim_{x_{out}\to\infty};}}[-x_{out},x_{out}]=(-\infty,+\infty), and (ii) the half-line \raisebox{2.15277pt}{\scalebox{0.8}{\displaystyle\lim_{x_{out}\to\infty};}}[x_{in},x_{out}]=[x_{in},\infty), with the density vanishing at both boundaries. In case (i) we find

[TABLE]

where the squared parenthesis denote the difference in the value of the function at the extremes viz., $[f(x)]_{{-x_{out}}}^{{x_{out}}}=f({x_{out}})-f({-x_{out}})$. With the exception of a few cases, such as the Bessel process with constant drift (discussed in §IV), the asymptotic behavior of the oscillating eigenfunctions satisfies a Helmholtz equation. By taking the limit $x\to\infty$ in the Schrödinger equation (29), and hence considering \raisebox{2.15277pt}{\scalebox{0.8}{\displaystyle\lim_{x\to\infty};}}\,\mathcal{V}_{s}(x)=\lambda_{cr}, the asymptotic behavior of the solution $\tilde{u}(k,x)$ satisfies

[TABLE]

so that the generic asymptotic behavior can be written as

[TABLE]

with amplitudes $\mathcal{A}(k)$ and $\mathcal{B}(k)$. At this stage we define one eigendifferential for each class of free solutions, whereas Fues defines eigendifferentials of linear combinations of solutions, which do not reproduce the solution of the Brownian motion problem discussed in §III, and additionally the associated integrals are very complicated.

Now, substituting the cosine in Eq. (49) into Eq. (47) the quantity in the square brackets in the integrand becomes

[TABLE]

where $w_{i}\equiv\sqrt{2\,k_{i}}/\sigma$, $i\in\{1,2\}$, and $\mathcal{A}_{i}\equiv\mathcal{A}(k_{i})$333Note that using the sine in Eq. (49) gives the same expression as the cosine, but leads to an opposite sign in the second row of Eq. (II.3), which eventually vanishes.. Using $\cos(a)\sin(b)=\frac{1}{2}[\sin(a+b)-\sin(a-b)]$, the previous equation becomes

[TABLE]

In order to integrate Eq. (47) over $k_{2}$ we note that the amplitudes $\mathcal{A}_{i}$ are very slowly varying functions of $k$, and thus we treat them as constants with respect to integration without introduction of significant error. The precision of this procedure must be thoroughly assessed numerically throughout the parameter space, making the method difficult to use.

Because Eq. (II.3) contains only odd functions so that $\sin(x)=-\sin(-x)$ and $\sin(0)=0$, we can apply the rule $[\,]_{-x_{out}}^{x_{out}}\to 2\,[\,]_{0}^{x_{out}}$, obtaining

[TABLE]

where

[TABLE]

where $\Delta_{n}k=[k_{in},k_{out}]$, and $\text{Si}(x)$ is the sine integral function. By taking the boundary to infinity, the cosines in Eq. (II.3) are suppressed while the sine integrals tend to a Heaviside theta function,

[TABLE]

where $\Theta(a+b)$ is zero (one) for $b<-a$ ($b>-a$). Due to the two Heaviside functions, for $k_{1}\in\Delta_{n}k$ and $k_{1}\in-\Delta_{n}k$, the integrals $I_{1}(k_{1},x_{out})$ and $I_{2}(k_{1},x_{out})$ are equal to $\pi\,\sqrt{k_{1}\,2}/\sigma$ and zero otherwise. Therefore, upon taking the infinite domain limit, the integral $I_{2}$ vanishes and the integral $I_{1}$ takes on values solely inside the shell, as depicted in Fig. 2. Therefore, the orthonormalization rule Eq. (41) is verified. Finally, we reparametrize the integral around $k_{1}=k_{2}$ as

[TABLE]

from which the squared normalization is

[TABLE]

The calculation for the half domain case proceeds analogously and, due to the different extremes in Eq. (II.3), differs only by a factor two;

[TABLE]

III Brownian motion with constant drift

In this section we demonstrate the technique by solving the Fokker-Planck equation for Brownian motion with a constant drift function $f(x)=a$

[TABLE]

and natural boundary conditions, $\rho(x=\pm\infty,t)=0$. Regardless of whether we have positive ($a>0$) or negative ($a<0$) deterministic drift, the potential $V(x)=-a\,x$ will be unbounded at $x=\pm\infty$. Therefore, the stationary solution

[TABLE]

is not normalizable, and the process is intrinsically transient. We can solve the Fokker-Planck equation by solving the associated Schrödinger eigenvalue problem, Eq. (29), with constant Schrödinger potential

[TABLE]

for which the Sturm-Lioville problem has only a continuum of eigenvalues. Since the stationary solution is not normalizable the eigenvalue $\lambda=0$ does not belong to the spectrum, which consists only of oscillating solutions.

We begin by writing the general solution $u_{\lambda}(x)$ of the Schrödinger equation,

[TABLE]

with $\text{C}_{1}$ and $\text{C}_{2}$ the integration constants. We exclude the case where $\lambda<a^{2}/(2\,\sigma^{2})$ because both solutions diverge at one boundary and vanish at the other. For $\lambda\geq a^{2}/(2\,\sigma^{2})$ both the linearly independent solutions are imaginary exponentials, so that sines and cosines are solutions viz.,

[TABLE]

with $\lambda=a^{2}/(2\,\sigma^{2})+k$. Note that these oscillatory eigenfunctions of Brownian motion have the same asymptotic behavior as Eq. (49), with amplitudes $\mathcal{A}=\mathcal{B}=1$. Therefore, the norm is given by Eq. (56) and the normalized eigenfunctions are

[TABLE]

The solution can then be written as a sum of two classes of eigendifferentials $\Psi_{n}^{cos}(x)$ and $\Psi_{n}^{sin}(x)$ as

[TABLE]

where

[TABLE]

and

[TABLE]

Here the shell size range is $\Delta_{n}k=[k_{n}-\eta,k_{n}+\eta]$ and the two sets of eigendifferentials are shown in Fig. 3 for different shell sizes.

Finally, the solution of the process can be determined by substituting Eq. (63) into Eq. (35), which gives

[TABLE]

By using the trigonometric relation

[TABLE]

and

[TABLE]

we obtain the well-known Gaussian distribution;

[TABLE]

with a broadening enhanced by the constant deterministic drift force characterized by $a$.

To determine the solution in the positive half domain, the general solution, Eq. (62), must satisfy the boundary condition at the origin. Hence, we discard the cosines, and the normalization is now provided by Eq. (57). However, because the remaining set of sines have a non-vanishing derivative at the origin, we have an absorbing boundary condition there and a natural boundary condition at infinity, giving the solution as

[TABLE]

which can be simplified using the integral relation

[TABLE]

finally leading to the solution of Eq. (58) as,

[TABLE]

We conclude this section by emphasizing that reproducing known results acts as the ideal test bed for the technique. Note that Eq. (73) can be obtained by employing the method of images and we also direct the reader to the work of Linetsky VALI2005 regarding the spectral representation of Brownian motion with a reflecting boundary condition at the origin (see also Abate and Whitt ABWH1987a ; ABWH1987b ) .

IV Bessel process with constant drift

Now we treat the Bessel process with constant negative drift, described a Fokker-Planck equation with potential $V(x)=-b\,\ln(x)-a\,x$, $a<0$ and $b>0$, with a natural boundary condition at infinity and a reflecting (i.e., no-flux) boundary at the origin that implies $\phi(x=0)=0$. The process has the normalizable stationary solution,

[TABLE]

where

[TABLE]

Thus, the stationary solution has a Bessel-like polynomial behavior near the origin, and is regularized at large $x$ by the exponentially-decaying tail associated with the linear term in the potential.

The Sturm-Liouville equation is closely related to the Schrödinger equation governing the radial component of the wave function of an electron with energy $E$ in an atom of atomic number $Z$, namely

[TABLE]

where $r$ is the radial coordinate, $y(r)=r\,R(r)$, with $R(r)$ the radial component of the wavefunction, and $l$ is the angular quantum number (c.f., Eq. 2.1 of Bethe and Salpeter (BSQM1957, )). The associated Schrödinger potential is

[TABLE]

which as seen in Fig. 4 tends asymptotically to a critical value

[TABLE]

defining the boundary between the discrete ($\lambda<\lambda_{cr}$) and the continuous ($\lambda>\lambda_{cr}$) spectra.

The Schrödinger potentials of the Bessel process of Eq. (4) and the Rayleigh process of Eq. (7), both with $b\to-b$, are

[TABLE]

As shown in Fig.5, whether the potentials and the Schrödinger potentials decay or diverge asymptotically leads to either continuum or discrete states. A peculiarity of the Bessel process is that, due to a scaling symmetry of the potential, its eigenfunctions are not oscillatory as we might expect, but instead decay at infinity ESGR2005 .

Now we describe the technique of solving the Schödinger-like equation for the Bessel process with constant drift. We note that when reparameterizing the equation in terms of Whittaker functions (see Appendix B), it describes an equivalent process with multiplicative noise SBMS1979 . The generic solution $u_{\lambda}(x)$ of the Schrödinger equation reads

[TABLE]

where $M$ and $W$ are the Whittaker functions of the first- and second-kind, $\text{C}_{1}$ and $\text{C}_{2}$ are the integration constants, and

[TABLE]

By using the relations

[TABLE]

and

[TABLE]

we can rewrite the solution as

[TABLE]

where $M$ and $U$ are the hypergeometric confluent functions of the first- and second-kind respectively, or the Kummer ($M$) and Tricomi ($U$) functions. Now, by using the symmetry property

[TABLE]

we can write the $W$-functions as

[TABLE]

IV.1 Discrete spectrum

By taking the limit $\lambda\to 0$ of Eqs. (IV) and (86), and considering that

[TABLE]

we find that only the $W$-function provides the stationary solution,

[TABLE]

where $U(0,\beta;x)=1$. Hence, from here forward we set $\text{C}_{1}=0$.

Upon imposition of the boundary conditions, $\psi_{\lambda}(x=0)=\psi_{\lambda}(x=\infty)=0$, we note that the condition at infinity is always fulfilled, whereas that at the origin is only satisfied when the first argument is a negative integer, thereby defining the discrete set of eigenvalues as

[TABLE]

with $n$ a non-negative integer so that

[TABLE]

For negative values of its first argument, the Tricomi function reduces to the unassociated Laguerre polynomial viz.,

[TABLE]

Substituting Eq. (90) into the expression (81) of $\zeta_{\lambda}$ gives

[TABLE]

and substituting this into $u^{(2)}_{\lambda}(x)$, Eq. (86), and using (91) leads to the expression

[TABLE]

in which

[TABLE]

The normalization $\text{C}_{2}$ is determined by the integral

[TABLE]

which gives

[TABLE]

Finally, the eigenfunctions in the discrete spectrum region are

[TABLE]

which form an infinite set (hence $N=\infty$ in the spectral representation of Eq. (II.2)) that are orthonormalized in the sense Eq. (33). However, the complete set is the sum of the eigenfunctions of the discrete and continuous spectrum444Completeness can also be inferred by reparameterizing the solution with $\tilde{x}=\xi_{n}\,x$ and $\gamma=-1+(2\,b)/\sigma^{2}$ obtaining $\psi_{n}(\tilde{x})=\mathcal{C}\,\omega(\tilde{x})\,L_{n}^{\gamma}(\tilde{x})$, with $\omega(\tilde{x})=\text{e}^{-\tilde{x}/2}\,\tilde{x}^{(\gamma+1)/2}$ and $\mathcal{C}$ a constant. The weight function of a complete basis of Laguerre polynomials is instead $\omega(\tilde{x})=\text{e}^{-\tilde{x}/2}\,\tilde{x}^{\gamma/2}$. In other words, a complete orthonormal basis of Laguerre polynomials is found by letting $\left(\xi_{n}\,x\right)^{\frac{b}{\sigma^{2}}}\to\left(\xi_{n}\,x\right)^{\frac{b}{\sigma^{2}}-\frac{1}{2}}$ in Eq. (97)..

IV.2 Continuous spectrum

In the continuous regime $\lambda>\lambda_{cr}=a^{2}/(2\,\sigma^{2})$ we are only concerned with the boundary condition at the origin; $\psi_{\lambda}^{con}(x=0)=0$. Note that for $\lambda>\lambda_{cr}$ the parameter $\zeta_{\lambda}$ defined in Eq. (81) takes on complex values. Therefore, we take only the real part of the complex eigenfunctions and parameterize the eigenvalues as

[TABLE]

with $k\geq 0$. To impose $\psi_{\lambda}^{con}(x=0)=0$ we study the behavior of the general solution at small $x$. Both solutions diverge at $x=0$ for $b>\sigma^{2}$ and vanish for $b<\sigma^{2}$. In the former case we have

[TABLE]

Thus, by taking a linear combination of the $W$- and $M$-functions we define a modified Whittaker function, $K$, in which the divergencies at the origin cancel each other;

[TABLE]

and thus $K$ vanishes at the origin, where $\alpha$, $\beta$ and $\gamma$ are the arguments of the functions in Eq. (99).

Finally, Eq. (57) gives the normalization of the eigenfunctions, which are

[TABLE]

where $\Re$ is the real part and $\mathcal{Y}(k)$ is the amplitude of asymptotic oscillations. Note that Eq. (101) also holds when $b<\sigma^{2}$. The evaluation and precision of the normalization coefficient are described in Appendix C.

IV.3 Numerical comparison

The spectral expansion of the transition density is

[TABLE]

where $\rho_{st}(x)$, $\psi_{n}^{dis}(x)$ and $\psi_{k}^{con}(x)$ are Eqs. (74), (97) and (101) respectively. Note that our expression is equivalent to that provided by Linetsky (see Eqs. 6 and 8 of LINE2004 ). However, while the form of the eigenfunctions of the discrete spectrum coincides with his, the continuum is written in terms of different Whittaker functions; complex $M$-functions in Linetsky’s work and a linear combination of $W$ and $M$-functions in our work. Of course the two formulations must be equivalent. In the left panel of Fig. 7 we compare our solution for the probability density to a numerical solution evaluated using the Crank-Nicholson method with a time step of $\Delta t=10^{-5}$ and an initial condition of $\rho(x_{0},t_{0})\propto x_{0}^{\beta}\,\text{e}^{\,\alpha\,x_{0}}$, where $\alpha=-16$ and $\beta=3$.

V Conclusion

The technique described in this paper is motivated by a stochastic model in climate dynamics, which treats the seasonal evolution of the Arctic sea ice thickness. The summer is modeled by a Brownian motion process with a negative constant drift and the winter by a Bessel process with a constant negative drift. We have solved the associated Fokker-Planck equations whose spectral expansions exhibit a continuum or a mixed (discrete and continuum) spectrum of eigenvalues. Relying upon early work in quantum mechanics, we posited that the normalization of oscillating eigenfunctions in the continuous spectrum depends solely upon the amplitude of the asymptotic oscillations, which led to a general formula for the solution of the one-dimensional autonomous Fokker-Planck equation with additive noise. We have employed this technique to solve the problem of Brownian motion with constant drift and natural boundary conditions at infinity and an absorbing boundary at the origin. Such well-trodden ground served as a good test bed for the approach, which we then demonstrated by solving the more complicated equation arising from the Bessel process with negative constant drift. Our solution coincides with that of Linetsky who employed different methods VALI2004 ; LINE2004 . However, the simplicity of our method of normalizing the functions in the continuum comes at the price of a lower computational efficiency. Therefore, depending on the specifics of the application, one may consider both methods on equal footing. Our approach can be used to examine the fully time dependent climatological evolution of the system by matching the two processes continuously, which is the subject of a future study.

Acknowledgements.

We thank S. Agarwal, F. Mancarella, M. Testa and S. Toppaladoddi for useful discussions. W.M. acknowledges a NASA Graduate Research Fellowship and a Herchel-Smith postdoctoral fellowship for support. W.M. acknowledges a Herchel-Smith postdoctoral fellowship for support. J.S.W. acknowledges NASA Grant NNH13ZDA001N-CRYO, Swedish Research Council grant no. 638-2013-9243, and a Royal Society Wolfson Research Merit Award for support. All authors acknowledge the 2015 Geophysical Fluid Dynamics Summer Study Program at the Woods Hole Oceanographic Institution, which is supported by the National Science Foundation and the Office of Naval Research under OCE-1332750.

Appendix A Determination of prefactor $\Delta_{n}k$

We start by considering the normalization relation for eigendifferentials,

[TABLE]

We then use the definition in Eq. (39) to rewrite the eigendifferentials in terms of eigenfunctions. Performing the spatial integral first, we obtain

[TABLE]

By using Eq. (34) we arrive at

[TABLE]

When evaluating the integral in $k_{2}$ for $m\neq n$, since $k_{1}$ lies outside of the range of integration of $k_{2}$, the delta function gives zero, thereby reproducing the behavior of the Kronecker delta. If instead $m=n$, the integration gives

[TABLE]

Appendix B Parameterization into the Whittaker equation

Here we parameterize the Sturm-Liouville equation with the potential in Eq. (77) into the Whittaker equation,

[TABLE]

which has the general solution

[TABLE]

where $\text{C}_{1}$ and $\text{C}_{2}$ are integration constants, and $M_{\{k,n\}}(x)$ and $W_{\{k,n\}}(x)$ are Whittaker functions of the first- and second-kind respectively, and are related to the confluent hypergeometric functions of the first- and second-kind via Eqs. (82) and (83). The confluent functions are related to the unassociated Laguerre polynomials through

[TABLE]

where $n$ is a positive natural number and $k$ and $x$ are real numbers. Finally we note that the $W$-function satisfies the identify $W_{\{\alpha,\beta\}}(x)=W_{\{\alpha,-\beta\}}(x)$.

We rescale Eq. (B.1) with $x\to\alpha\,x$, so that the equation and its solution are now

[TABLE]

Multiplying both sides with $\beta$, setting $\alpha=\sqrt{2\,\beta}/\sigma$ and rescaling $k\to k/\sqrt{\beta}$ we find

[TABLE]

Finally, we obtain Eq. (77) by letting

[TABLE]

Appendix C Normalization of the modified

Whittaker function, $K$.

Here we evaluate the asymptotic amplitude of the oscillating solutions in the continuous spectrum, and numerically test the accuracy of Eq. (57). We start by taking the asymptotic behavior of Eq. (100) in a polar representation viz.,

[TABLE]

where

[TABLE]

with

[TABLE]

in which $\Gamma$ is the Euler Gamma function, $\csc$ the cosecant, and $\Im$ and $\Re$ are the imaginary and real parts respectively. We can sum the four terms as

[TABLE]

where

[TABLE]

and

[TABLE]

in which $\text{arctan}_{2}$ is the two-argument arctangent. Therefore, the asymptotic behavior of the solution is

[TABLE]

and the coefficient in Eq. (57) is the normalization only if $\phi_{1}$ varies with $k$ more rapidly than does the amplitude $\mathcal{Y}$ and the other phases. As shown in Fig. 8, the approximation is robust, save for the small region near $k=0$, that contributes negligibly to the full integral over $k$.