{'fraction_non_alphanumeric': 0.044900760739003606, 'fraction_numerical': 0.012326259795229652, 'mean_word_length': 4.447421716778314, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solution of the Fokker-Planck equation with a logarithmic potential and mixed eigenvalue spectrum

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.


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. 1,2]. Statistical information is governed by the probability density function ρ(x, t), which satisfies the Fokker-Planck equation a second order parabolic partial differential equation wherein f (x) describes the deterministic dynamics and σ is the constant noise strength. (Note, for parsimony of notation, we write the operator ∂ 2 /∂x 2 as ∂ 2 x .) For example, when a one-dimensional solely diffusive process governs ρ(x, t), there is no deterministic drift force (f (x) = 0), and the Fokker-Planck equation is a diffusion equation with a Gaussian solution; ρ(x, t) = 1 (2 π σ 2 t) for initial data ρ(x, t = 0) = δ(x), with natural boundary conditions ρ(x = 0, t) = ρ(x = ∞, t) = 0 1 . The long-time behavior of Eq. (1) for f (x) = 0 is given by the Boltzmann distribution in the case of thermal systems, defined as where N is a normalization constant and we can consider V (x) = − 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 which has the following exact solution in which I 1 2 + b σ 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,∞), the origin satisfies a reflecting boundary condition and there is a natural boundary condition ρ(x = ∞, t) = 0 at infinity.
The Bessel process exhibits several interesting properties. For example, as described in [refs. 4, 5], by integrating over angle variables and letting b → b−(d−1) σ 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 ind = d − 2 b/σ 2 dimensions, whered can be negative and non-integer. Pitman and Yor [6,7] 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., 8,9]. Hence, squared Bessel processes play a major role in the study of maps between Brownian motion in different dimensions [10,11].
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 [12]; (b) the vortex-antivortex annihilation in the 2D XY Heisenberg model subject to thermal fluctuations at temperatures below the Kosterlitz-Thouless transition T KT [4]; (c) an alternative technique to quantize non-Abelian gauge field theories, avoiding the drawbacks of the Fadeev-Popov method [13] 2 . 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 [14,15].
As summarized nicely by Dechant et al., [16], 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 ρ st (x) ∝ e − 2 σ 2 V (x) = x −2 b/σ 2 clearly diverges at large x (small x) depending on whether 2 b/σ 2 < 1 (2 b/σ 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) ≈ 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 [17], the potential may also exhibit a more complex logarithmic dependence, such as 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 ρ st , but instead requires taking the infinite-time limit of the time-dependent moments [16]. The latter has been proved by Lutz to be related to a breaking in the ergodicity for systems with power-law distributions [18]. By adding regular terms to the potential, the existence of the stationary solution can always be guaranteed. One example is the Rayleigh process [19][20][21][22], given by 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 [23]. The Rayleigh process is widely used in finance because the parabolic term corresponds to the Cox-Ingersoll-Ross (CIR) process [24]. 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 [15].
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, 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 [25]. 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 [25]. 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. [26] 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 [27]. Linetsky [28] 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 [29,30]. The concept of the sea ice thickness distribution was introduced by Thorndike et al. [31], 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 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 ρ(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 [32] 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 [29,33,34]. 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 [35], 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 δ-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 [27], and by solving the equation in complex space in the full domain [0, ∞) [28] (see [36] 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 where 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 Ω = [x in , x out ]. We will consider the simpler case of additive noise, where g(x) ≡ σ 2 defines the constant noise strength, and time-independent constant drift f (x, t) ≡ f (x), and hence the Fokker-Planck equation becomes corresponding to which is the autonomous Langevin equation where dW (t) is a Wiener process. The probability density can be written as where ρ 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 lim t→t0 ρ(x, t|x 0 , t 0 ) = δ(x − x 0 ), and thus formally Thus, as described in §I, the stationary solution of Eq. (11) is where N is a normalization constant, viz., and we can consider V (x) = − x f (y) dy to be a potential associated with the deterministic drift force; . 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.

A. 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 [37]. We expand the transition density as where the φ λ are the set of eigenfunctions of L satisfying the eigenvalue problem and Λ is the eigenvalue spectrum. Inserting the expansion (17) into the density (13) we find where the coefficients α λ (t 0 ) are the projections of the initial condition on the λ th function, viz., 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 L in Eq. (18) is not self-adjoint, which can be seen by writing L in terms of the potential V (x) as and checking that for two test functions φ 1 and φ 2 the following inequality holds Nevertheless, we can construct a self-adjoint operator associated with the Fokker-Planck operator as which satisfies Eq. (18) and possesses a complete orthonormal set of eigenfunctions {ψ λ } with non-negative real eigenvalues λ ≥ 0. The two operators have the same spectrum, and their eigenfunctions are, up to normalization, related by The completeness of the basis insures that the initial condition of the transition density satisfies which, by adding unity, in the sense of distributions, can be rewritten as Substituting (26) into (14), and then into (19) using (24), we find B. 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 which resembles an imaginary-time Schrödinger equation with potential Again, the Sturm-Liouville theorem guarantees that the set of eigenfunctions, ψ λ (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 Λ d , (b) become a continuum Λ c , or (c) be a mix of the two Λ = Λ d + Λ 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 operatorL is self-adjoint, the spectrum in this region can be solved simply by imposing the boundary conditions on the generic solutions ψ λ (x) = C u λ (x), where C is an integration constant. The discrete set of solutions satisfying the boundary conditions will occur at values λ n , n = 0, 1, 2, · · · , with λ 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 = n, namely, where n ≡ λ n and m ≡ λ m are two eigenvalues. The integration constant is determined by normalization of the solution as to satisfy the orthonormality relation where δ 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, ψ 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 δ-function as 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 Finally, at one (or both) of the boundaries approaching infinity, the potential may asymptotically approach a critical value λ 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 ψ dis λn , with eigenvalues λ n , n = 0, 1, 2, · · · for λ < λ cr , and a continuous set, ψ con k , with eigenvalues k ∈ [0, ∞), where λ = λ cr + k. Thus, solution for the transition density in this case is where N is the highest discrete eigenvalue for a finite discrete set.

C. Normalization of oscillating functions
As mentioned in §II B, 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 → ∞ as in [38,39]. Here we adapt the approach introduced in the 1920s by Fues [40,41] 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 in which N and λ cr are the same as in Eq.(36), a k (t 0 ) is the projection of the initial condition on the k-th eigenfunction, and φ(x) ∝ ρ st (x) ψ(x). We then split the integral into a series of infinitesimal integrals viz., where the ∆ n k, n = 1, 2, · · · , ∞, form an infinite discrete set of eigenvalue shells. Here, k n denotes the eigenvalue of the n th shell such that lim ∆nk→0 k n = k, and the functions Ψ n (x) are the eigendifferentials, defined as the integrals of eigenfunctions over the n th shell; Although a single eigenfunction ψ con k (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, for all n and m. Thus, we can generalize the orthonormality definition to Determination of the factor ∆ n k in the orthonormalization rule for Ψ n (x) is given in Appendix A. The projection a k (t 0 ) in Eq. (37) can then be defined as Now, by using the eigendifferential definition of Eq. (39) in the third of Eqs. (41), we obtain for n = m where {k 1 , k 2 } ∈ ∆ n k. The prefactor (∆ n k) −1 insures that the above integral holds for any arbitrary shell size so long as where u(k, x) is the general solution of the Sturm-Liouville equation and 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 ψ con k (x). Let us take the Schrödinger-like equation (29) for two different eigenvalues k 1 and k 2 , multiply them by ψ con k2 and ψ con k1 respectively, and then subtract them, which gives Now, using Eq. (45) in Eq. (44), and dividing by (k 1 − k 2 ), leads to We integrate by parts in a finite domain with the boundary terms defined at infinity, and consider two domains; (i) where the squared parenthesis denote the difference in the value of the function at the extremes viz., [f (x)] xout −xout = 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 → ∞ in the Schrödinger equation (29), and hence considering lim x→∞ V s (x) = λ cr , the asymptotic behavior of the solutionũ(k, x) satisfies so that the generic asymptotic behavior can be written as with amplitudes A(k) and 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 where w i ≡ √ 2 k i /σ, i ∈ {1, 2}, and A i ≡ A(k i ) 3 . Using cos(a) sin(b) = 1 2 [sin(a+b)−sin(a−b)], the previous equation becomes In order to integrate Eq. (47) over k 2 we note that the amplitudes 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. (51) contains only odd functions so that sin(x) = − sin(−x) and sin(0) = 0, we can apply the rule where and where ∆ n k = [k in , k out ], and Si(x) is the sine integral function. By taking the boundary to infinity, the cosines in Eq.
(53) are suppressed while the sine integrals tend to a Heaviside theta function, . Due to the two Heaviside functions, for k 1 ∈ ∆ n k and k 1 ∈ −∆ n k, the integrals I 1 (k 1 , x out ) and I 2 (k 1 , x out ) are equal to π √ k 1 2/σ 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 from which the squared normalization is The calculation for the half domain case proceeds analogously and, due to the different extremes in Eq. (51), differs only by a factor two;

In this section we demonstrate the technique by solving the Fokker-Planck equation for Brownian motion with a constant drift function
and natural boundary conditions, ρ(x = ±∞, 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 = ±∞. Therefore, the stationary solution 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 for which the Sturm-Lioville problem has only a continuum of eigenvalues. Since the stationary solution is not normalizable the eigenvalue λ = 0 does not belong to the spectrum, which consists only of oscillating solutions. We begin by writing the general solution u λ (x) of the Schrödinger equation, with C 1 and C 2 the integration constants. We exclude the case where λ < a 2 /(2 σ 2 ) because both solutions diverge at one boundary and vanish at the other. For λ ≥ a 2 /(2 σ 2 ) both the linearly independent solutions are imaginary exponentials, so that sines and cosines are solutions viz., with λ = a 2 /(2 σ 2 ) + k. Note that these oscillatory eigenfunctions of Brownian motion have the same asymptotic behavior as Eq. (49), with amplitudes A = B = 1. Therefore, the norm is given by Eq. (56) and the normalized eigenfunctions are The solution can then be written as a sum of two classes of eigendifferentials Ψ cos n (x) and Ψ sin n (x) as where Ψ cos and Here the shell size range is ∆ n k = [k n − η, k n + η] and the two sets of eigendifferentials are shown in Fig. 3 for different shell sizes.
Ψn(x) Finally, the solution of the process can be determined by substituting Eq. (63) into Eq. (35), which gives By using the trigonometric relation we obtain the well-known Gaussian distribution; 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 which can be simplified using the integral relation finally leading to the solution of Eq. (58) as, 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 [42] regarding the spectral representation of Brownian motion with a reflecting boundary condition at the origin (see also Abate and Whitt [43,44]) .

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 φ(x = 0) = 0. The process has the normalizable stationary solution, where 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 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 [41]). The associated Schrödinger potential is which as seen in Fig. 4 tends asymptotically to a critical value defining the boundary between the discrete (λ < λ cr ) and the continuous (λ > λ cr ) spectra.
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 [45]. 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 [46]. The generic solution u λ (x) of the Schrödinger equation reads where M and W are the Whittaker functions of the first-and second-kind, C 1 and C 2 are the integration constants, and By using the relations and we can rewrite the solution as 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 we can write the W -functions as A. Discrete spectrum By taking the limit λ → 0 of Eqs. (84) and (86), and considering that we find that only the W -function provides the stationary solution, where U (0, β; x) = 1. Hence, from here forward we set C 1 = 0. Upon imposition of the boundary conditions, ψ λ (x = 0) = ψ λ (x = ∞) = 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 with n a non-negative integer so that λ n = (a n σ) 2 + 2 a 2 b n 2 (b + n σ 2 ) 2 .
For negative values of its first argument, the Tricomi function reduces to the unassociated Laguerre polynomial viz., Substituting Eq. (90) into the expression (81) of ζ λ gives and substituting this into u (2) λ (x), Eq. (86), and using (91) leads to the expression in which .
The normalization C 2 is determined by the integral which gives Finally, the eigenfunctions in the discrete spectrum region are which form an infinite set (hence N = ∞ in the spectral representation of Eq. (36)) that are orthonormalized in the sense Eq. (33). However, the complete set is the sum of the eigenfunctions of the discrete and continuous spectrum 4 .

B. Continuous spectrum
In the continuous regime λ > λ cr = a 2 /(2 σ 2 ) we are only concerned with the boundary condition at the origin; ψ con λ (x = 0) = 0. Note that for λ > λ cr the parameter ζ λ defined in Eq. (81) takes on complex values. Therefore, we take only the real part of the complex eigenfunctions and parameterize the eigenvalues as with k ≥ 0. To impose ψ con λ (x = 0) = 0 we study the behavior of the general solution at small x. Both solutions diverge at x = 0 for b > σ 2 and vanish for b < σ 2 . In the former case we have 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; and thus K vanishes at the origin, where α, β and γ are the arguments of the functions in Eq. (99). Finally, Eq. (57) gives the normalization of the eigenfunctions, which are where is the real part and Y(k) is the amplitude of asymptotic oscillations. Note that Eq. (101) also holds when b < σ 2 . The evaluation and precision of the normalization coefficient are described in Appendix C.

C. Numerical comparison
The spectral expansion of the transition density is where ρ st (x), ψ dis n (x) and ψ con k (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 [28]). 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 ∆t = 10 −5 and an initial condition of ρ(x 0 , t 0 ) ∝ x β 0 e α x0 , where α = −16 and β = 3. The plot was evaluated using 20 eigenfunctions in the discrete spectrum and an AccuracyGoal equal to 3 for the continuum. Because our solution is written in terms of a complicated asymptotic linear combination, Mathematica took approximately five times longer in absolute time to evaluate the integral over k relative to the solution of Linetsky [28].

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 [27,28]. 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.
in which arctan 2 is the two-argument arctangent. Therefore, the asymptotic behavior of the solution is lim x→∞ [K] = Y(k) cos(φ 1 (k, x) + φ 2 (k, x) + φ 3 (k) + W(k)) (C. 6) and the coefficient in Eq. (57) is the normalization only if φ 1 varies with k more rapidly than does the amplitude 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.