Brownian motion in trapping enclosures: steep potential wells, bistable wells and false bistability of induced Feynman–Kac (well) potentials

We investigate signatures of convergence for a sequence of diffusion processes on a line, in conservative force fields stemming from superharmonic potentials U(x) ∼ x m , m = 2n ⩾ 2. This is paralleled by a transformation of each mth diffusion generator L = DΔ + b(x)∇, and likewise the related Fokker–Planck operator L* = DΔ − ∇[b(x) ⋅], into the affiliated Schrödinger one Ĥ=−DΔ+V(x) . Upon a proper adjustment of operator domains, the dynamics is set by semigroups exp(tL), exp(tL*) and exp(−tĤ), with t ⩾ 0. The Feynman–Kac integral kernel of exp(−tĤ) is the major building block of the relaxation process transition probability density, from which L and L* actually follow. The spectral ‘closeness’ of the pertinent Ĥ and the Neumann Laplacian −ΔN in the interval is analyzed for m even and large. As a byproduct of the discussion, we give a detailed description of an analogous affinity, in terms of the m-family of operators Ĥ with a priori chosen V(x)∼xm , when Ĥ becomes spectrally ‘close’ to the Dirichlet Laplacian −ΔD for large m. For completness, a somewhat puzzling issue of the absence of negative eigenvalues for Ĥ with a bistable-looking potential }0, m{ >}2$?> V(x)=ax2m−2−bxm−2,a,b,>0,m>2 has been addressed.


Motivation
Our present investigation is motivated by that of reference [1], where a number of puzzling features of Lévy jump-type processes in bounded domains, in particular these with reflecting boundaries, were reported. The concept of reflected Lévy flights has found a coverage in the mathematical literature [2]. However, if the fractional dynamics (e.g. Lévy-stable or fractional Brownian motion) is to be restricted to a bounded domain, in view of the nonlocality of the process and its generator, there are many inequivalent proposals for fractional analogs of the domain-restricted Laplacian (Dirichlet, Neumann and the likes).
This ambiguity is particularly annoying in the reflective case. It shows up in attempts to invent an appropriate process behavior at the boundary and extends to its particular technical implementation, i.e. not only to identifying a proper nonlocal version of the Neumann condition, but also to a detailed path-wise definition of the reflection proper at the boundary. Examples are provided by computer-assisted procedures for the generation of sample paths and varied realizations of the reflection or quasi-reflection (happening in a vicinity of the boundary) scenario [3][4][5]. The inferred jump-type processes are inequivalent as well, differing in their spectral and statistical characteristics. This affects theoretical predictions for relaxation properties and the near-equilibrium asymptotic behavior.
On the other hand, while departing from the Langevin picture of the Lévy motion that is confined by superharmonic potentials ∼ x m /m, m = 2n, n 1, one arrives at fractional Fokker-Planck equations, whose limiting properties as m → ∞ (the convergence issue is mathematically disputable), are ultimately interpreted in terms of reflected Lévy flights in the interval, [6][7][8]. The point is that no link has been established with the Neumann fractional Laplacian (neither with any convincing form of the Neumann condition), which is a valid generator of the Lévy process in a bounded domain with reflecting boundaries.
In our opinion, the pertinent link can be consistently analyzed by employing the transformation of the fractional Fokker-Planck equation to the fractional Schrödinger-type equation. It is a properly tailored version of the technical tool, often used in the study case of the standard (Brownian) Fokker-Planck equation, but seldom addressed in the literature on confined Lévy processes. Since the superharmonic approximation seems to provide a suggestive method to understand what is possibly meant by the reflected Lévy process in the bounded domain (we restrict considerations to the interval on R), the pertinent transformation to the Schrödingertype dynamics should in principle provide an approximation of this governed by the Neumann Laplacian.
We emphasize that our intro on conditioned Lévy processes stands merely for a departure point and actually provides motivation for the present Brownian endeavor. From now on, we shall not refer to the Lévy framework anymore. Instead, we turn back to the extremally confined Brownian motion, where the sequential limiting behavior (superharmonic approximation, with a reflecting Brownian motion introduced as a formal limit) is taken for granted, or in the least unproblematic.
A transformation of the Fokker-Planck operator to the Schrödinger one (a respected metod of solution for the Fokker-Planck equation, [25,26]) might seem to be 'obvious' as well. The expected outcome of the emergent sequential approximation is a possibly 'smooth' limiting behavior of Schrödinger-type operators −Δ + V(m) as m 1 and their spectral data, while approaching these normally attributed to the reflected Brownian motion, where the motion generator is the standard Neumann Laplacian (−Δ) N with the well known spectral solution [26][27][28][29][30][31][32][33][34].
Links with the semigroup picture (Schrödinger-type dynamics) are technically much simpler to analyze in case of the confined Brownian motion, than in case of Lévy flights. Therefore, we decided to re-address this issue in some detail, with the aim of detecting any hitherto overlooked/neglected subtleties and obstacles/deficiencies of the formalism (including computerassisted procedures) for: (i) superharmonic confining conditions, (ii) permanent confinement between reflecting boundaries, (iii) mutual relationships between (i) and (ii), with an emphasis on the issue of the 'spectral closeness' of the superharmonic Brownian motion generator to the standard Laplacian with Neumann boundary data and the reliability of involved approximations. And, likewise on the 'closeness' (whatever that means, but with the path-wise picture in mind) of the superharmonically confined Brownian motion to the standard reflected one.
One-dimensional diffusion processes, which are confined in the interior of the finite interval in R, are typically classified to be reflecting (Neumann boundaries) or taboo (Dirichlet boundaries which are not accessible), and directly involve spectral solutions of Laplacian operators restricted to that interval, and respecting the boundary data, i.e. Δ D and Δ N respectively [25,26,30,31].
For each boundary data choice, the lowest (interval-restricted) Laplacian eigenvalues (strictly speaking, the spectral gap) determine the time-rate of the so-called spectral relaxation to an invariant pdf (probability density function) of the Fokker-Planck equation in question. In the Dirichlet case, the drifted diffusion in the interval derives from the killed one via a proper conditioning [26,30,32,33], and the bottom eigenfunction and the lowest eigenvalue of Δ D are explicit constituents of the analytic formula for the transition probability density of the process with inaccessible boundaries.
In the literature there exist many proposals on how to approximate Brownian motion in the interval by means of strongly confining random model systems, both on the Langevin (thus Fokker-Planck) and Schrödinger type (semigroups, generalized diffusion equations) levels of description. Typically one employs either a sequence of finite well potentials with an increasing height of the barrier, or a sequence of strongly anharmonic potentials, alternatively ∼ x m , x m /m, m x m with m = 2n growing indefinitely (plus a number of other options) [34][35][36][37][38][39][40][41].
These procedures tell us how to get at a reliable approximation of the Dirichlet (infinite well/finite interval-restricted) operator (−Δ) D , in terms of Schrödinger-type generators and semigroup evolutions (with the Feynman-Kac paths redistribution involved [17]), [30][31][32][33][34][35][36][37][38] see also [39]. Although there still remain mathematical obstacles that annoy physicists [40,41]. In this case, a parallel sequential approximation in terms of affiliated Langevin-Fokker-Planck evolutions is unavailable in the literature and this topic, together with a number of illuminating explicit examples, we shall address in the appendices.
In comparison with the Dirichlet case, the Neumann one has been left somewhat aside. In the current literature, it is taken for granted (see [1] and references there-in) that the Brownian motion in extremally anharmonic force fields (Langevin-driven case with b(x) ∼ x 2n−1 , n 1) not only provides a satisfactory approximation of the reflected Brownian motion in the interval, but actually converges to the latter.
At this point we indicate a problem that will appear in our subsequent analysis, that the convergence is not unifom, and that the standard Neumann condition is actually not reproduced in the controlled point-wise limit. This statement needs to be reconciled with the fact that the analytic description of the reflected motion in the interval is given in terms of the spectral solution for the Laplacian with Neumann boundary conditions, Δ N , see [27][28][29].
We point out well established mathematical facts [32]: (i) the Schrödinger operator with a constant potential and the Neumann boundary condition corresponds to the reflected Brownian motion, while (ii) the Schrödinger operator with a constant potential and the Dirichlet boudary condition corresponds to what we name the taboo process [30], e.g. a drifted Brownian motion that is conditioned never to exit the open interval in question. Moreover, (iii) there is a one-to one correspondence between the pertinent Schrödinger operators (thence semigroup dynamics) and drifted diffusion generators (thence Langevin-driven diffusion process).
This motivates the basic aim of the present paper, to re-examine the theory of the Brownian motion on the interval, in terms of a sequence of steep (superharmonic) potential wells, with an emphasis on the familiar and widely exploited [25,26] link between the Fokker-Planck dynamics ∂ρ = L * ρ and that governed by the affiliated Schrödinger semigroup exp(−tĤ ), e.g. the generalized diffusion equation ∂ t Ψ = −ĤΨ. The latter equation is known to quantify the random motion in terms of the Feynman-Kac formula [52][53][54][55][56], and might be interpreted path-wise as the Brownian motion 'in potential landscapes', see [17], see also [22][23][24] for an alternative interpretation of the 'landscape' notion.
Our departure point will be the case of Brownian motion in conservative force fields (explicitly present in Langevin and Fokker-Planck equations), which stem from confining potentials of the form U(x) = x m /m, m = 2n, n 1. Asymptotic stationary probability densities (invariant pdfs of Fokker-Planck equations), to which the random dynamical system relaxes, have the Boltzmann form ρ * (x) ∼ exp(−U(x)/D).
Although the local minima of our bistable potentials have values that are below zero, negative energy bound states are conspicuously absent. At the first glance, this property seems to contradict standard intuitions underlying the concept of bistability and the related tunellingthrough-the-barrier conceptual imagery of quantum theory (cf a 'canonical' discussion of the negative eigenvalues splitting in case of the familiar double-well potential [50]). This slightly puzzling point will be addressed in below as well.

Boltzmann-type equilibria for Smoluchowski diffusion processes
Let us consider a one-dimensional diffusion process [25], with the Langevin representatioṅ where B(s) = 0, B(s)B(s ) = δ(s − s ) and b(x) is a forward drift of the process having the gradient form b = 2D∇Φ, D being a diffusion constant. If an initial probability density ρ 0 (x) is given, then the diffusion process obeys the Smoluchowski-Fokker-Planck equation where L * = DΔ − ∇(b·) is the Fokker-Planck operator (we prefer the notation L * instead of the often employed L FP ). We recall [25,26], that L * is a Hermitian L 2 (R) adjoint of the diffusion operator L = DΔ + b∇, which in the mathematical literature actually stands for a legitimate generator of the pertinent Markovian stochastic process. Let us introduce an osmotic velocity field u = D ln ρ and the current velocity field v = b − u. The latter directly enters the continuity equation ∂ t ρ = −∇j, where j = v · ρ has a standard interpretation of a probability current.
We restrict further discussion to time-independent drifts, that are induced by external (conservative, Newtonian) force fields f = −∇U. One arrives at Smoluchowski diffusion processes by setting This expression accounts for a fully-fledged phase-space derivation of the spatial (Smoluchowski) process, in the large friction β regime. It is taken for granted that the fluctuationdissipation balance gives rise to the standard form D = k B T/mβ of the diffusion coefficient [25,35].
In the stationary asymptotic regime we have j → j * = 0. We denote ρ * = ρ * (x) a strictly positive probability density, to which ρ(x, t) is presumed to relax as t → ∞. Accordingly in Consequently, we have a stationary solution of the Smoluchowski (Fokker-Planck) equation where 1/Z is a normalization constant.
Here F * = −k B T ln Z is the minimal value of the time dependent Helmholtz free energy of the random motion F = F(t) = U + k B T ln ρ , to which F(t) relaxes as t → ∞ [35], . denotes the mean value with respect to ρ(x, t).
Note that we consider confining potentials U(x) with the property lim |x|→∞ U(x) = +∞ and such that exp(−U(x)/k B T ) is L 1 (R) integrable (to secure the normalization). Any bounded from below potential function U(x) can be redefined to become either non-negative or positive, since an additive energy 'renormalization' is irrelevant for the Langevin-Fokker-Planck dynamics and leaves the the functional form of ρ * (x) intact (this in view of the L 1 (R) normalization).

Pseudo-Schrödinger reformulation of the Fokker-Planck dynamics
For computational convenience it is useful to scale away the factor mβ in the formulas of the previous subsection. This can be done by changing the time scale in the Fokker-Planck equation (2) from t to t/mβ, accordingly replace D by D ≡ k B T, and thus b( We shall intentionally employ the notation D ≡ k B T, instead of D = mβD = k B T. From now on our D in reality will be a 'graphical' replacement for D (with time rescaling kept in memory).
Accordingly, if we consider literally equation (2)  Although in principle D > 0 is unrestricted, in below we shall mostly refer to the moderate noise range and specifically to D = 1 or D = 1/2 (that to stay in conformity with mathematically oriented approaches to diffusion problems). This, to some extent, may blur a direct comparison with observable effects of the specific noise intensity choice. Indeed, small D narrows/compresses the potential profile, while increasing its steepness. The choice of larger D flattens/stretches the potential profile). An issue of the small noise D 1, and small versus moderate or large noise, will briefly return back in section 4.
The outlined procedure may be regarded as a reconstruction of the semigroup dynamics from an eigenstate (here, ground state function) [19]. It is an alternative to another, reverse engineering procedure [21,42] (named also a targeted stochasticity), whose main goal is to reconstruct the random motion (e.g. the Langevin-driven motion and the Fokker-Planck equation) that is compatible with a priori given equilibrium function (e.g. the stationary pdf) and a priori prescribed noise (here Brownian/Wiener one).

Brownian motion in potential landscapes (Feynman-Kac potential spatial profiles)
All previous observations actually derive from the Langevin-type equation (1) and thus from the traditional (microscopically motivated) picture of the Brownian motion in conservative force fields. Forces are of external nature, while Wiener noise is regarded as an intrinsic property of the environment in which the motion takes place.
Concerning the generalised diffusion equation equation (7), let us set D = 1/2 and make an assumption that the potential function V is a continuous, bounded from below function. Then, we can introduce the positive symmetric integral kernel k(t, x, y) = k(t, y, x) of the semigroup operator exp(−tĤ ), given e.g. by the Feynman-Kac formula, with an explicit V entry, see [30,55,56], see also [1,30].
The Feynman-Kac formula defines the transition density (integral kernel of the semigroup operator) as a weighted integral over sample paths of the diffusion process (with the Wiener path measure being involved; here paths ω originate from x at time 0 and their destination is y to be reached at time t > 0): In reference [17] equation (11) is interpreted in terms of a random mover in an energy landscape set on R by the potential V(x) and acting as a mechanism that reinforces or penalizes the random mover tendency to reside or go into specific regions of space. A 'responsibility' for a weighted redistribution of random paths in a given time interval, is here transferred from external force fields of equation (1) to the spatial variability of potentials V(x) of equation (8), in particular their curvature and steepness.
The Fokker-Planck-induced dynamics of ρ(x, t), equation (2), in view of equations (6) and (7) can be interpreted as follows: which in a self-defining manner introduces [30], the transition probability density p(t, y, x) of the (conditioned) Markovian diffusion process underlying the temporal evolution of ρ(x, t).
We emphasise the relevance of the order in which variables x and y appear, because p(t, x, y) = p(t, y, x). Indeed, in contrast to the kernel function k(t, x, y), transition pdfs are not symmetric functions of x and y, which has profound consequences in the derivation of the explicit functional form for generators L and L * of the diffusion process, while evaluated according to the familiar recipes. Namely the semigroup T t = exp(tL) has a generator The time evolution of probability measures and associated probability density functions ρ(x, t) is governed by the adjoint semigroup While coming back to the unrestricted noise intensity parameter D, we realize that the pdf dynamics (2) and (11), as a direct consequence of equation (7), can be represented in various forms [20,[22][23][24]: Remembering that ρ * (x) ∼ exp[−U(x)/D], we uncover a direct impact of the Boltzmann weight upon the dynamics of ρ(x, t). It suffices to set ρ ±1/2 * → exp[∓U(x)/2D] everywhere in equation (11).
In the Langevin modeling of the Brownian motion, the choice of a Newtonian potential U(x), determines both an invariant pdf ρ * = (1/Z )exp(−U/D) of the F-P equation and the forward drift (e.g. the conservative force field) b(x) = D∇ ln ρ * = −D∇U, entering ∂ t ρ = DΔρ − ∇(bρ). The relaxation process ρ(x, 0) → ρ(x, t) → ρ * (x) in question is well posed, once a suitable ρ(x, 0) is selected as the initial data.
In the dual picture provided by the Feynman-Kac (Schrödinger) semigroup, an a priori chosen ρ * (x) determines the Feynman-Kac potential V = D[Δρ 1/2 * ]/ρ 1/2 * and thence the 'potential landscape' delineated by the spatial profile of V [17]. The relaxation process refers to the time rate at which the eigenfunction ρ Remark: In the study of an impact of strongly inhomogeneous (disordered) media upon the random motion (with an independent noise input) [24], a complementary notion of salience (attractivity) field has been associated with the exponential (Boltzmann) weight . A curve delineated by the salience field s(x) (actually its square root in equation (12)), might deserve as well an interpretation (modulo dimensional issues) of the 'energy (salience) landscape'. It is the relative weight between the source and target locations that enforces or penalizes (reduces) the rates of spatial displacements of the random process.

Joint spectral solution for motion generators L, L * andĤ
As long as the operators L, L * ,Ĥ are considered in a common for them function space L 2 (R), the transformation from L * to the Hermitian (symmetric, eventually self-adjoint) operatorĤ seems to be the natural way to address spectral properties of generators of motion and related relaxation time rates to equilibrium (we need the real eigenvalues ofĤ ) at the bottom of the non-negative spectrum). It is however not a must and we can admit other function spaces.
Given a unique stationary solution (5) of the Fokker-Planck equation, L * ρ * (x) = 0. The uniqueness is granted if the Fokker-Planck operator has an energy gap at the bottom of its spectrum [26,32,33]. This gap in turn controls the (generically exponential) time rate of relaxation of ρ(x, t) toward the equilibrium pdf ρ * (x). This is the scenario of so-called spectral relaxation.
Concerning the existence of the gap, a sufficient condition to this end can be given by demanding that the Gibbs potential U(x) is a continuous and continuously differentiable (e.g. C 2 (R)) function, and that the semigroup potential V(x), equation (8) shares with U(x) the confining property: lim |x|→∞ V(x) = +∞, see e.g. [26].
Independently of the previous factorization (6) of ρ(x, t), given ρ * (x), we can introduce an alternative factorization: A comparison with equation (6) shows that Ψ(x, t) = h(x, t)ρ 1/2 * (x). A direct consequence of the Fokker-Planck equation (often named the forward Kolmogorov equation) is that h(x, t) is a solution of the backward Kolmogorov equation: where L = DΔ + b∇ stands for the diffusion generator and the initial condition reads h( It is not yet widely appreciated information that the generator L is Hermitian (symmetric [32], and eventually self-adjoint) in the Hilbert space L 2 (R, ρ * ), whose scalar product is weighted by the invariant pdf ρ * : [26]. We point out that −L is a nonnegative operator, with a discrete spectrum beginning form the bottom eigenvalue 0.
We can readily introduce an operator product (warning: we admit action upon functions in the domain in a definite order, from right to left), reproducing the operatorĤ of equation (11): We point out that the spectrum ofĤ, −L and −L * is nonnegative. This is accompanied by the transformations: and L = ρ −1 * L * ρ * . Three operatorsĤ, L and L * are Hermitian (and eventually self-adjoint) in function spaces L 2 (R), L 2 (R, ρ * ) and L 2 (R, ρ −1 * ) respectively [26]. The transformations (14)-(16) allow to deduce eigenfunctions and eigenvalues of each related spectral problem, if one of them has been actually solved.
Namely, let us assume that we know the spectral data ofĤ. We denote ψĤ k the kth eigenfunction in L 2 (R) and λĤ k the associated eigenvalue. This eigenvalue (up to a sign, (15)- (17)) is shared by all three generators, while the corresponding eigenfunctions are related accordingly: (19) and In particular we have: We recall that the pertinent eigenfunctions 'live' in different function spaces and any reduction to L 2 (R) (unweighted) integrals needs to be executed with due caution.
The above transformations reduce the issue of the spectral analysis of the Brownian motion to that of the computational convenience. See e.g. [26-29, 32, 33] for eigenfunction expansions associated with the diffusion generator L of a reflected process on the interval.
The main value of the spectral analysis of motion generators is that we may compare their spectral gaps, if in existence, under various physically motivated circumstances. They allow to quantify the convergence to equilibrium for an ergodic diffusion process and actually coincide with the exponential time rates of convergence.

Spectral versus nonspectral relaxation in the OU process. A comment
At this point it is worthwhile to invoke a discussion of a possibility of the non-spectral relaxation in the one-dimensional Ornstein-Uhlenbeck process [43]. The basic argument therein can be rephrased as follows (we employ our notation). Equation (6) is supposed to define a similarity transformation between the Fokker-Planck operator L * andĤ, cf equations (16)- (18) (remember about the operator action ordering from the right to the left).
If we have the spectral solution forĤ in hands, in terms of eigenvalues λ and L 2 (R) eigenfunctions Ψ λ (x), then equation (17) tells us that the eigenvalues of L * are −λ, while the eigenfunctions appear in the transformed form φ λ (x) = Ψ λ (x)ρ 1/2 * (x). The probability distribution ρ(x, t) that can be expanded into φ λ (x) will relax at rates, determined by gaps in the eigenvalue spectrum ofĤ. This is called a spectral relaxation pattern.
It is clear that the validity of the aforementioned expansion is guaranteed only, if L * is a Hermitian operator. To this end an appropriate function space needs to be specified, and we know that ρ(x, t) (and likewise the initial datum ρ(x, 0)) must be an element of L 2 (R, ρ −1 * ), as mentioned before in a paragraph below equation (8). This excludes from considerations such initial pdf as the Cauchy distribution, mentioned in reference [43] (the pertinent discussion appears between the formulas (3) and (4)), in connection with the conjectured non-spectral relaxation pattern.
The point is that to employ the Fourier series (eigenfunction) expansion one needs to respect carefully adjusted domain properties of transformed operators, see (15)- (20). Our freedom of choice of initial conditions for ρ(x, t) is thus limited to these pdf ρ(x, 0) only, which decay sufficiently faster at infinity than ρ −1/2 * (x) grows, and thus transform via equation (6) into square integrable functions Ψ(x, t) [26]. Told otherwise, the Cauchy pdf is not a legitimate candidate for an initial distribution for the Fokker-Planck equation of the Ornstein-Uhlenbeck process.
Assume that ρ * (x) stands for the invariant pdf of the Ornstein-Uhlenbeck process (cf the appendices section for more details). Let f n (x) stand for the eigenfunctions of the OU process generator L = Δ − x∇ in L 2 (ρ * ), with eigenvalues −λ n = n, n = 0, 1, 2, . . . . Then, according to equation (20), L * (ρ * f n ) = −nρ * f n . In terms of the eigenfunctions of L * , we have the the expansion formula: Consequently, beginning from an arbitrary initial distribution ρ 0 (x) ∈ L 2 (R, ρ −1 * ), the probability law of the process converges exponentially fast to the invariant distribution. This is the spectral relaxation scenario and there is no room for its violation in the theoretical setting of the standard Ornstein-Uhlenbeck process..

Extremally steep (superharmonic) potentials in the Langevin-Fokker-Planck approach
It is a folk wisdom, that a sequence of symmetric single well (superharmonic) potentials U 2n (x) = (x/L) 2n /2n, with the growth of their steepness (i.e. for large values of n), can be used as an approximation of the infinite well potential with reflecting boundaries located at x = ±L. The statement has been widely adopted in the literature, its validity taken for granted and possible limitations have been ignored or bypassed. While some caution is here necessary.
For computational simplicity, in below, we shall pass to the dimensionless notation x/L → x, so that U(x) = x 2n /2n, for large values of n stands for an approximation of the infinite well, supported on [−1, 1]. Closely related potentials of the form U(x) = 2n x 2n and U(x) = x 2n likewise can be employed as the infinite well approximations.
We point out that one needs to observe possibly annoying boundary subtleties. Namely, at x = ±1 the pertinent potentials take values 1/2n, 2n and 1 respectively, and their respective limiting values (point-wise limits, as 2n grows to ∞) read 0, ∞ and 1. That needs to be kept under scrutiny, while invoking the standard definition of the infinite well enclosure, set on the interval [−1, 1] ⊂ R (we impose the exterior boundary data, instead of locally defined data which are normally assigned to the interval endpoints only), see [39,40]: This infinite well definition stays in conformity with the limiting n → ∞ behavior of U(x) = 2n x 2n at x = ±1, but is incompatible with the remaining U(x) options, which take boundary values 1 or 0, instead of ∞.
Recently, the 2n → ∞ sequential approximation has been explored in the Langevin-Fokker-Planck description of the Brownian (and subsequently Lévy-stable) motion in the extremally anharmonic (Newtonian) force field b(x) ∼ −x 2n−1 , to mimic and ultimately justify the validity of reflecting boundary data in the emergent infinite well enclosure [13][14][15][16]. The affiliated Schrödinger-type route has never been analyzed.
On the other hand, an analogous limiting procedure, while employed on the Schrödingertype operator level of description, has been found to provide an accurate approximation of the traditional Dirichlet infinite well enclosure [36][37][38][39]. Surprisingly, in this case the Langevin-Fokker-Planck alternative has not been explored at all (this point we shall address in section 5). Let where the integral expression is recognizable as the Euler Gamma function . Accordingly, we have: The normalization coefficient B for the case of U m (x) = x m , m = 2n, is the m 1/m multiple of that in equation (22): B = m 1/m A. Since m 1/m → 1 with m → ∞, the limiting behavior (convergence rate) in both cases is similar.
We note that A ≡ A m approaches the limiting value 1/2 from below as a growing function, while B ≡ B m approaches 1/2 from above as a decreasing function.  Since ρ(0) = A or B respectively, we realize that the pertinent normalization coefficients actually set the value of an asymptotic 'plateau', with 1/2 referring to a uniform probability distribution on the interval [−1, 1], formally regarded as the m → ∞ limit for the sequential approximation x m /m. However, the uniform distribution is recovered on an open interval (−1, 1) only, if one would employ x m and mx m approximating sequences. See e.g. figures 1-3.
It seems instructive to visualize the dependence of ρ * (x) on m = 2n in the close vicinity of boundary points x = ±1. This is depicted in figure 4 for U(x) = x m /m, at points x = 0.99, 1, 1.01.

The reference problem: reflected Brownian motion on [−1, 1]
The uniform probability distribution on the interval [−1, 1] may be interpreted as a signature of reflecting (Neumann) boundary data and thence of the reflecting Brownian motion in the interval [26][27][28]. We emphasize that nothing is normally said about the exterior R\[−1, 1] of the pertinent interval and in the traditional mathematical lore, the Neumann data are defined locally [26,32].
An exact solution of the latter problem refers to the standard Laplacian, while restricted to the interval and subject to reflecting (e.g. Neumann) boundary conditions. Solutions of the diffusion equation The pertinent transition density reads [28]: and is an integral kernel of the reflecting semigroup exp(tΔ N ). The Neumann operator −Δ N admits the eigenvalue 0 at the bottom of its spectrum, the corresponding eigenfunction being a constant 1/ √ 2, whose square actually stands for a uniform probability distribution on the interval of length 2, see [26,28,32].
Solutions of the diffusion equation with reflection at the boundaries of D = [−1, 1] can be modeled by setting p(x, t) = k N (t, x, x 0 ), while remembering that p(x, 0) = δ(x − x 0 ). We can as well resort to Ψ(x, t)= D k N (t, x, y)Ψ(y) dy. Note that all n 1 eigenvalues coincide with these attributed to the absorbing case [30], and (up to dimensional constants) coincide with eigenvalues (nπ/2) 2 of the standard (quantum mechanical) infinite well problem, with Dirichlet boundary data. The eigenfunctions respect Neumann conditions, and do not necessarily vanish at the boundary points (that would be the Dirichlet case).
We note that the exponential rate of convergence to equilibrium in the present case (D = 1) reads π 2 /4 and is three times smaller than the rate E 2 − E 1 = 3π 2 /4 established in section A.4 as an asymptotic signature of the taboo process (conditioned never to exit the interval with inaccessible boundaries [30,32]).
The direct path-wise description of the reflected Brownian motion belongs to a non-standard inventory, if compared with the standard Langevin modeling. It involves the so-called Skorokhod problem and a class of stochastic differential equations with reflection [28,29]. This problem is avoided (or circumvented) in the pragmatic approach to the reflection issue via computer simulations of sample paths, where it is the boundary behavior (proper handling of the instantaneous reflection) which appropriately alters the statistical features of propagation of the otherwise free Brownian motion [4,6,15].

Superharmonic approximations of the reflecting well: (in)validity of the sharp Neumann boundary condition
As long as we prefer to deal with traditional Langevin-type methods of analysis, it is of some pragmatic interest to know, how reliable is an approximation of the reflected Brownian motion in [−1, 1] by means of the attractive Langevin driving (and thence the Fokker-Planck equation), with force terms (e.g. drifts) coming from extremally anharmonic (steep) potential wells.
The main obstacle, we encounter here is that a 'naive' m = 2n → ∞ limit is singular and cannot be safely executed on the level of diffusions proper. We note that for any finite m, irrespective of how large m actually is, we deal with a continuous and infinitely differentiable potential and likewise, the Boltzmann-Gibbs pdf as a consequence of (1)-(5).
On the informal, graphical level (figures 2 and 3) we can anticipate that the limiting pdf, is a constant 1/2 in the closed interval (uniform distribution on the interval) and vanishes identically for |x| > 1: This is consistent with the formal limiting behavior of U(x) = U m (x) = x m /m (with U m explicitly present in the exponent of the BG density), as m → ∞, and properly mimics the infinite well potential enclosure [8], whose boundary ±1 can be continuously reached from the interval interior: We note that the infinite-valuedness everywhere beyond [−1, 1] (e.g. an exterior restriction) is an essential complement to potential shapes depicted in figures 2 and 3.
If (25) is taken literally as a (more or less mathematically legal) limit of a sequence of superharmonic potentials, there is an obstacle to be overcome or bypassed. Namely, differentiability properties of resultant (limiting) functions U and ρ * are lost at the interval [−1, 1] endpoints. We emphasize that an infinite differentiability is a valid property of U(x), ρ * (x) and ρ 1/2 * (x) for all 2 m < ∞, irrespective of how large m actually is.
Let us consider in more detail the properties of −∇U and ∇ρ * in R, for large m = 2n, in the vicinity of ±1. Our focus is on the seemingly obvious reflecting behavior and the possible (in)validity of Neumann boundary conditions, which actually define the reflecting Brownian motion in terms of the Neumann Laplacian Δ N in the interval, cf subsection 4.2.
In the superharmonic Langevin-type regime we have b(x) = −∇U(x) = ∇ ln ρ * (x) = −x m−1 . Since ∇ρ * (x) = −Ax m−1 exp(−x m /m), we readily infer the location of its minimum in the vicinity of x = 1: which, for large m, can be safely replaced by x = m 1/m . One may infer, see [57], a useful estimate for the actual location of the considered minimum for finite but large values of m: with a limiting property lim m 1/m = 1 as m → ∞.    Accordingly for all finite values of m, irrespective of how large m is, there is a lot to happen (in the lore of turned over sample paths) in the narrow zone of thickness 4/ √ m, beyond the interval boundaries set at ±1 on R.

Reconstruction of the Schrödinger semigroup from an eigenstate of the motion generator
While inspired by the targeted stochasticity concept of reference [42], we follow a procedure described in [20] to infer the Schrödinger semigroup from a given a priori invariant probability density function (strictly speaking, from its square root). The semigroup dynamics is to deduced from the knowledge of the strictly positive zero energy eigenfunction ρ /m, m = 2n 2 of the sought for (semigroup) motion generator, see figure 7 and references [19,20,42] for rudiments of the reverse engineering concept and this of the targeted stochasticity.
By setting D = 1 in equations (7)- (10), we interpret the introduced m-family of ρ withĤ = −Δ + V, where the reconstructed (given ρ 1/2 * (x)) Feynman-Kac potential reads, cf figures 8 and 9: Remark: We recall that the effects of noise (with an arbitrary intensity D > 0) can be reintroduced through the formula m) . A specific choice of the the noise intensity value becomes immaterial in the large m limit.
In the literature, slightly less restrictive reconstruction problems for the Schrödinger-type dynamics have been considered. Namely [19,44], one may depart as well from the eigenvalue problem in L 2 (R) with the given a priori L 2 (R) eigenfunction ψ 0 (x) (typically free of nodes, albeit it is not a must, see remark in below), and ask for the well behaved potential function which makes the eigenvalue problem forĤ properly defined and solvable. This is in fact the main idea behind the 'reconstruction of the dynamics form the eigenstate' [19,54]. It is clearly an inversion of the standard logic, where one first selects the appropriate potential and then seeks solutions (and eigensolutions) of the Schrödinger or Schródinger-type (generalized diffusion one, our case) equation. Remark: It is customary to reconstruct the dynamics from the strictly positive function, which is interpreted as the ground state of the sought for motion generator. However, it is not a must. One may as well admit other eigenfunctions, that lead to non-negative probability densities, instead of the strictly positive one. It is possible to handle the problem of zeroes. A detailed analysis of the harmonic oscillator case in this regard (evaluation of drifts and related semigroup potentials of the form (8) from excited eigenstates) is provided in the text between formulas (46) and (50) of reference [52]. See also reference [53].
The functional form (43) of V(x) looks quite intriguing, if compared with the familiar anharmonic expression |x| m , notoriously employed in the quantum theory literature as the potential function for the m-family of Schrödinger operators, whose spectral properties in the m → ∞ limit have been a subject of elaborate studies [38].
In the present paper, we are mostly interested in the large m behavior of the potential (43), and specifically whether the reconstructed semigroup spectral problem might be perceived as 'close' to that related to the Neumann Laplacian Δ N , compare e.g. section 4.2.

On the spectral affinity (m 2 regime) with the reflecting Brownian motion on the interval
Let us consider the eigenvalue problem where V(x) has the two-well functional form (30). We are interested in testing whether with the growth of m one may relate the corresponding semigroup generator with the Neumann Laplacian Δ N in L 2 ([−1, 1]). Let us choose the Neumann basis on [−1, 1], see [1,30] and section 4.2: ψ n (x) = cos nπ 2 (x + 1) , n = 0, 1, . . .
Presuming that it is the reference basis system, we make a hypothesis that any eigenfunction ψ(x) in equation (32), for large m, should be 'close' to the corresponding Neumann eigenfunction. Likewise, we expect the same 'closeness' property for the corresponding eigenvalue (in the large m regime). If (32) is a valid eigenvalue problem, the eigenvalues should follow from the L 2 (R) expectation value λ n,m = ψ n |Ĥ|ψ n ,  that needs to be approximated by the expectation value ofĤ in L 2 ([−1, 1]) with respect to the Neumann basis.
According to figure 8, for |x| > 0.99 the singular behavior of V(x), in the large m regime, is incompatible with the Neumann basis properties at the interval boundaries. If the boundary parameter a goes to 1, we need to take into account the term (m − 1)x m−2 , which does not vanish at ±1 and gives a large input at the boundary, in the numerically assisted integration procedure.
Therefore, instead of extending the involved integration to the whole of [−1, 1], we restrict the interval boundaries to [−a, a], a < 1. Accordingly, we disregard the misbehaving part of the integral (this misbehavior is a consequence of replacig/apprixmating the 'true' eigenfunction by the Neumann one). Moreover, we isolate a leading term (nπ/2) 2 , n 1 which is characteristic for the Neumann well spectral solution. The remaining integral expression is expected to decay to 0 as m → ∞:  We can give another argument supporting the statement that signatures of the aforementioned spectral affinity in the interval [−a, a], a 0.99 have been met, by depicting ( figure 13) the detuning diagram, i.e. the plot of |Hψ 1 − λ 1,m ψ 1 | for a = 0.9 and m = 10, 30, 100, as a function of x, provided ψ 1 (x) ≡ cos(π(x + 1)/2) in [−a, a]. The detuning can be made arbitrarily small with the growth of m, as long as we stay within [−a, a], a < 1.

False bistability of two-well Feynman-Kac potentials
Rather specific functional form of the deduced semigroup (Feynman-Kac) potential (30) and its increasingly singular behavior at the boundaries of [−1, 1], with the growth of m, makes problematic a direct computer-assisted evaluation of higher level eigenvalues and eigenfunctions. The level of involved technical difficulties has been noted before [44], in connection with an example of the sextic potential V(x) = ax 6 − bx 2 . We note that for m = 4, e.g. U(x) = x 4 /4, (43) reduces to the sextic potential with a = 1/4 and b = 3/2. For m = 6, e.g. U(x) = x 6 /6, (43) takes the form of the decatic potential V(x) = ax 10 − bx 4 with a = 1/4 and b = 5/2.
We mention an extended discussion of decatic potentials presented in references [44,46,47], which can be regarded as a reference record of various curiosities concerning the solvability of spectral problems closely related to ourĤ = −Δ + V(x). These curiosities are shared by a broad class of quasi-exactly solvable Schrödinger equations, whose solutions are not amenable to standard algebraic methods and basically unknown, except for very few examples. A direct recourse to computer assisted arguments is unavoidable.
A related problem is the false bistability of bistable-looking two-well Feynman-Kac potentials. Namely, as follows from the previous discussion, operatorsĤ = −Δ + V(x), with V defined by equation (30), have a nonnegative energy spectrum beginning from the eigenvalue zero. As noticed in reference [44], in this case traditional WKB approximation methods are inapplicable to retrieve spectral data, because the classical zero energy level in the considered case refers to an unstable equilibrium point of the two-well potential.
Moreover, an insight into computer-generated solutions shows that the system appears to ignore ('does not feel') the presence of two potential wells, contrary to expectations rooted in the folk understanding of diffusion problems in double-well potentials [25,26]. (Note that we make a distinction between two-well and double-well potentials; the latter have always a quartic form).
Let us mention the classic double-well potential: U(x) = ax 4 /4 − bx 2 /2, with a > 0, b > 0. The potential has two negative-valued minima (local well depth) at x 1,2 = ± b/a, and a local maximum equal 0 at x 0 = 0. (Note that a modified potential U(x) + b 2 /4a = a(x 2 /2 − b/2a) 2 is nonnegative and has a local maximum b 2 /4a, while local minima take the value zero). In the course of the Langevin-driven Brownian motion with a drift b(x) ∼ −bx 3 + ax, even a weak noise causes particles not to remain in the same (local) well, and makes them to 'oscillate' forth and back between two wells. That should possibly result in the bimodal stationary pdf, with two conspicuous maxima and a minimum in between them (possibly at x 0 ).
One might justifiably be tempted to invoke the standard imagery the escape time analysis (for a randomly wandering particle from the well of temporary residence etc.) [58][59][60][61][62][63]. However, as follows from our further discussion (section 5.2) this temptation is somewhat deceiving and relies on the relative balance between parameters a and b.

When do we encounter the bottom eigenvalue zero for two-well systems, or why there is no signatures of the potential bistability?
Double-well potentials in the classic quartic version have received an ample coverage in the literature, specifically in connection with effects of tunneling through the barrier, which gives rise to a splitting of otherwise degenerate negative energy (might be positive, but then necessarily with energies below the top of the central double-well barrier i.e. below the local maximum of the potential). A standard procedure, employed to analyze the low-lying spectrum, amounts to a local aproximation of each well by a suitable harmonic oscillator potential [50].
Let us consider the general eigenvalue problem for a Hermitian (eventually self-adjoint) Schrödinger type operator −Δψ + Vψ = λψ. (36) We presume the potential to be bounded from below and symmetric (this assumption does not harm the generality of further arguments). Let V(x) has a minimum (in the least a local minimum) at x 0 ∈ R and consider a Taylor expansion of V(x), which in a sufficiently close vicinity of x 0 may be reduced to Let us furthermore assume that V(x 0 ) < 0 in the vicinity of x 0 and at x 0 . Inserting an approximate expression (32) for V(x) to the eigenvalue equation (33), we arrive at an approximate eigenvalue problem: This expression can be readily compared with the standard harmonic oscillator spectral problem. Indeed, since we have: Equations (52) and (53) While keeping in mind that 2 2m ≡ 1, we get For the standard quantum harmonic oscillator the ground state eigenvalue n = 0 equals E 0 = ω/2, we can infer a corresponding λ 0 from the (approximate) identity E 0 = λ 0 − V(x 0 ), so arriving at An assumption that in the vicinity of the local minimum V(x 0 ), the potential function V(x) is negative, implies a rough existence criterion for negative eigenvalues in the two-well potential case. Indeed, a demand λ 0 < 0 enforces Both sides of the inequality are positive, hence after taking the square and employing (42) we can rewrite the existence criterion for the negative λ 0 , as the restriction upon the curvature of the two well potential V(x) at x 0 : If this inequality holds true, the existence of negative eigenvalues is allowed in the spectral solution forĤ. The above argument is not limited to the double-well (quartic) potentials and remains valid for more general two-well potentials. We work with a special subclass in the two-parameter family of two-well potentials: for which two local (well) minima are symmetrically located at ±x 0 , where Note that assuming b > 3a, for all m we have x 0 > 1, while for b < 2a the minima are located in the interior of [−1, 1]. Our specific case of the two-well potential equation (30) is identified by setting a = 1/4, b = (m − 1)/2 and assuming m = 2n, n > 1 (we recall that the case n = 1 refers to the standard harmonic oscillator potential with ground state energy subtracted [30,56]). This implies Further calculation is performed for our two-well potential (30). We have: and This implies a sharp inequality Since we assume m > 2, the condition m < 4/3 is violated from the start, which provides a formal explanation to why negative eigenvalues do not exist in the spectrum ofĤ = −Δ + V with the two-well potential V in the form (30). This conclusion stays in an obvious conformity with the fact that (according to the derivations of section A.2)Ĥ has the eigenvalue zero corresponding to the unimodal strictly positive eigenfunction. In turn, the latter property validates our harmonic oscillator approximation at local minima of the looking-bistable two-well potential and the inferred, equation (45) existence criterion for negative eigenvalues.   potentials, it may serve as a useful credibility test of our methods. Arguments of the previous subsection can be readily adopted to this familiar example. Let us consider and investigate how the the existence or non-existence of negative eigenvalues for H = −Δ + V depends on the mutual balance of steering parameters a and b. We note that the local maximum is located at 0, two negative-valued local minima (depth of the wells) equal −b 2 /4a at x 1,2 = ± b/2a, and the distance between the minima (e.g. an effective width of the central barrier) equals 2 b/2a. Note that a nonnegative definite version of the double-well potential In conformity with our previous discussion, let us address the existence issue of negative eigenvalues for the operatorĤ = −Δ + V, with V(x) given in the general double-well form (51). We employ the inequality (45), whose validity means that negative eigenvalues are permitted. We need to evaluate the local minimum value of V(x 0 ) and the curvature of the potential {V} (x 0 ). These read: V(x 0 ) = −b 2 /4a and {V} (x 0 ) = 4b. Accordingly, the negative eigenvalue existence criterion takes the form Let us examine this criterion for a simple but popular in the literature double well example: with two minima at ±α, α > 0 (the effective width of the cetral barrier equals 2α). The subtracted term α 4 actually sets the depth of each well (alternatively -height of the central barrier whose local maximum at x = 0 equals zero). A passage from the previous notation involves substitutions: a = 1, b = 2α 2 .
Accordingly, the bistability of two-well potentials (and the double-well in this number) may not be reflected in the topological properties of ground state eigenfunctions. These functions may be unimodal or bimodal, depending on the value of α, and there is a transitional value α critical which results in the eigenvalue zero of the operator −Δ + V.
For α < α critical the ground state eigenvalues are positive, while for α > α critical negative ground state eigenvalues are admitted. Thus, the topology change of the ground state function form unimodal to bimodal shape is reflected in the sign change of related eigenvalues ( figure 14).
Zero energy eigenvalue corresponds to the 'transitional shape', where the associated α critical stands for a bifurcation point: the local maximum of the unimodal function degenerates and bifurcates into twin local maxima of the bimodal one. The local maxima existence is the signature of bistability and the presence of negative eigenvalues forĤ = −Δ + V.

Reverse engineering for the Feynman-Kac quartic double well: reconstruction of the
associated Langevin-driven Brownian motion. It is figure A1 which, as a byproduct of the discussion of the double-well V(x), equation (51), shape impact upon the functional form (e.g. shape variability from modal to bimodal) of the ground state function and the bottom eigenvalue ofĤ, is our departure point. We have in hands all numerical data necessary to retrieve the affiliated Langevin-driven Brownian motion. We effectively follow the so-called reverse engineering procedure [20,21,42].
For each value of α, we employ numerical data specifying the related ground state function ρ 1/2 * (x). Accordingly, we have ρ * (x) in hands. This in turn allows to deduce (numerically as well) the corresponding drift field b(x) = −∇ ln ρ * (x) of the affiliated Langevin-Fokker-Planck driven dynamics. Given b(x), the functional form of nonegative valued (an additive constant input is here irrelevant) potentials U(x) have been reconstructed.
Accordingly, Langevin-driven Brownian motions with drifts b(x) depicted in figure B1, follow a spectral relaxation pattern to an invariant pdf ρ * (x), with time rate determined by the bottom eigenvalue of hatH which is shared with actual L and L * , see e.g. section 2.  figure A1). The Newtonian (Boltzmann, U = −ln ρ * ) potential U(x), such that b = −∇U, is depicted as well.

Specifying V(x) if U(x) is quartic.
We have mentioned before that the Langevindriven process with b(x) = −x 3 , i.e. U(x) = x 4 /4, is intimately related with the Brownian motion 'in a potential landscape' set by the Feynman-Kac potential in the sextic form V(x) = x 6 /4 − 3x 2 /2. It is interesting to know that a thorough dicussion of the mapping of the the Fokker-Planck evolution, with b(x) inferred from U(x) = ax 4 − bx 2 , to the corresponding pseudo-Schrödinger one (incorporating effects of arbitrary noise intensity D, has been accomplished in reference [62]. Namely, for U(x) = ax 4 − bx 2 , a, b > 0 and D > 0, we get the sextic Feynman-Kac potential: where: We note that the mutual balance of parameters a and b is reflected in the uni-or bimodality of the Fokker-Planck invariant pdf ρ * (x). Accounting for the presence of noise (e.g. accepting any D > 0) implies that the associated sextic Feynman-Kac potential may topologically vary from two-well to the threewell shape (this to be compared with an intuitive reasoning of reference [58]).

Sextic two-well potential.
In contrast to the quartic double-well, the sextic two-well potential upon identifications a = 1/4 and b = 3/2, becomes a member of our family (30). Coming back to the notation (55), we realize that the extrema are located at 0 and ±x 0 where 1/2 and V (x 0 ) = 8b. Accordingly, a criterion for the existence of negative eigenvalues takes the form 27a < b 2 .
In passing we note that in the case which corresponds to the bottom eigenvalue zero, the pertinent inequality does not hold true.
In reference [44] a one-parameter family of sextic potentials has been introduced (α > 0 is presumed). By inspection one can verify that for each member of this family the operator −Δ + V has the energy zero eigensolution of the form ψ 0 (x) ∼ exp(−αx 4 /4). We can readily verify that (identify a = α 2 and b = 3α) the sufficient condition for the existence of negative eigenvalues 27a < b 2 does not hold true. Plugging α = 1 we get V(x) = x 6 − 3x 2 with the related ψ 0 (x) ∼ exp(−x 4 /4). Note that the exponent x 4 /4 is twice the exponent x 4 /8 (inferred from ρ 1/2 * ), which we have associated with the potential function (30) In reference [44], and example is given of the sextic potential V(x) = x 6 − 7x 2 , for which by inspection we can verify that −Δ + V has two explicitly known eigenfunctions and eigenvalues. Namely: the ground state eigenfunction (yet unnormalized) ψ 0 ( 2) exp(−x 2 /4) corresponds to the positive eigenvalue E 0 = +2 √ 2 (it is not the first, but the second excited eigenfunction, has the same parity as the ground state). Plugging a = 1 and b = 7 we readily verify that the condition 27a < b 2 in the present case holds true, as expected.

Outlook
The transformation of the Fokker-Planck dynamics to the affiliated Schrödinger-type evolution has been analyzed in detail for a sequence of superharmonically confined systems, with a special emphasis on spectral properties of the Schrödinger semigroup and the role of its (Feynman-Kac) integral kernel in the construction of stochastic models of relaxation, stemming from the knowledge of the ground state function of the Schrödinger Hamiltonian and the integral (Feynman-Kac) kernel of the related Schrödinger semigroup.
The expected outcome should have been a possibly 'smooth' limiting behavior (this proved not to be the case, pointwise-limits do not reproduce the standard 'sharp' Neumann condition) of Schrödinger-type operatorsĤ = −Δ + V(m) as m → ∞, that would justify a reliable approximation of the reflected Brownian motion with the generator (−Δ) N by the Brownian motion in extremally steep anharmonic potentials. This (reliable approximation) property is normally expected to be transferable to the path-wise behavior as well.
In view of rather unusual two-well shape of pertinent Feynman-Kac potentials, and their singular at the boundaries m → ∞ limiting behavior, there is no smooth passage between two formalisms (e.g. an approximate and sharp model of reflection). Nonetheless, the approxima-tion by superharmonic models can be justified, if we admit a disregard of the narrow zone of a finite 'thickness' (cf our considerations in sections 3 and 4, and especially tables and figures therein), in the vicinity of the boundary. The optimal 'thickness' depends on how large m is.
We find worth emphasizing an analysis of section 5, where the concept of false bistability has been introduced. It is really amusing, that in the statistical physics literature nobody has ever paid attention to the generic issue of the eigenvalue zero of the diffusion generator (which, together with the Fokker-Planck operator is transformed into the affiliated HamiltonianĤ ), and non-existence of negative eigenvalues ofĤ for a broad family of two-well potentials, and the classic quartic double-well problem in particular. In this regard, we have deduced a nonexistence criterion, which extends earlier perturbative studies of the (quantum) double-well, [45].
A special role in the paper is played by appendices A-C. To enhance the coherence of our arguments, we have moved to appendix A a description of all consecutive steps being executed in the construction of transition pdfs of the conditioned Brownian motion (killing, its taming, and recovering Brownian relaxation for the Ornstein-Uhlenbeck process and the likes). Motion in the interval with Dirichlet boundaries is described in detail as well. In appendix B, we demonstrate the workings of specific sequential approximations, whose limiting case is the Brownian motion in the interval with inaccessible Dirichlet boundaries. Appendix C is of utmost importance, because there we demonstrate how the relaxation properties of conditioned Brownian motions can be deduced within the classic, but seldom revived in the literature, Schrödinger's boundary data and interpolation problem. That includes a derivation 'for nothing' of the basic functional form (C11) and (C12), shared by various transition probability densities introduced in the present paper, and derived by other methods (albeit with the Schrödinger problem reference) in reference [30].
The present Brownian endeavor has been planned as a detour, which actually turned over to an extended preparatory step, in the attempt to tackle theoretical and computer-assisted problems arising in the approximate description of Lévy-type processes in bounded domains with reflection. Partial results in existence, were obtained along Brownian-inspired lines (Langevin-Lévy motion in steep anharmonic potentials versus fractional Schródinger spectral problems) [64]. In passing, we note that the mapping of the fractional Fokker-Planck dynamics into the affiliated fractional Schródinger-type problem, has been investigated in the past [17][18][19][20][21][22][23][24], in connection with 'Lévy motion in energy landscapes'.
Leaving aside an approximation issue, reflected Lévy processes proper are nowadays an active research topic, albeit physical and mathematical papers on the subject show quite different, often incongruent, attitudes and proposals toward the proper description of (i) the involved path-wise behavior and (ii) a proper formulation of Neumann boundary condition for nonlocal operators [1,64].
Once more, we indicate that in the present paper we have addressed problems that are casually considered as 'obvious' and taken for granted, or possibly irrelevant in the considered contexts (this refers to the validity of the m → ∞ limit, the related approximation accuracy control, and to the nonexistence of negative eigenvalues for a broad class of two-well potentials, not forgetting about a profound problem of the quasi-exact solvability [44] of related Hamiltonian systems). The Brownian route has been chosen as a playground for checking jeopardies and possible inadequacies of the "obvious" strategy, and to 'clean the stage' prior to passing to an analysis of technically more involved issue of reflected Lévy flights and their superharmonic approximations, along the lines indicated in reference [64].

A.1. Ornstein-Uhlenbeck process vs harmonic oscillator semigroup: dimensional issues
Since, the explicit presence of dimensional constants somewhat blurs a connection between the Langevin equation-induced Fokker-Planck) dynamics (1)-(5), the inferred semigroup one (6)-(9) and the emergent spectral problem (10), we shall discuss in some detail the dimensional issues in the case of the standard harmonic attraction U(x) = kx 2 /2.
The drift b(x) = −(k/mβ)x = −κx defines the Ornstein-Uhlenbeck process, for which the semigroup (Feynman-Kac) potential (8) takes the form: A multiplication by 2mD restores the standard dimensional version (joule as the energy unit) of the potential, so we realize that upon setting ω = κ, the potential can be interpreted as that of a harmonic oscillator system, with an mDω subtracted. One should not confuse this ω with a natural frequency k/m associated with the harmonic potential U(x) = kx 2 /2. By means of a formal identification D ≡ /2m where is the reduced Planck constant, we give 2mDĤ the familiar form: of the quantum harmonic oscillator Hamiltonian with a ground state energy renormalization. More general discussion of this subtraction issue can be found in references [19,30,54,56], see also section 3 in reference [9] for spectral aspects of the relationship between the OU process and the harmonic oscillator semigroups. We know that the eigenvalue problem (9), while extended to the operator (A3), has the spectral solution in L 2 (R) with eigenvalues (2mD)λ n = E n − E 0 = ωn, n 0. The bottom eigenvalue equals zero. The corresponding L 2 (R)-normalized (ground state) eigenfunction reads: where U(x) = kx 2 /2, and to recover the functional form (5) of the invariant density ρ * (x), we have reintroduced the 'thermal' notation, e.g. ≡ 2mD, D = k B T/mβ and ω = κ = k/mβ. Accordingly F * = (1/2)ln(k/2πk B T ).

A.2. Eigenfunction expansions: killing can be tamed
We note that if the energy is measured in units of ω, while the distance in units of /mω, the rescaled energy operator with subtraction takes the form with the spectrum E n = n, n 0 and the ground state function φ 0 (x) = π −1/4 exp(−x 2 /2). Note that the spectrum of 2Ĥ = −Δ + x 2 − 1 coincides with 2E n = 2n and begins from the eigenvalue zero. Consider the (rescaled, without subtraction) harmonic oscillator problem denotedĤ 0 = (1/2)(−Δ + x 2 ). Its spectral solution comprises a sequence of eigenvalues n = n + 1 2 and corresponding eigenfunctions φ n (x) = [4 n (n!) 2 π] −1/4 exp(−x 2 /2) H n (x) which are L 2 (R) normalized. Here H n (x) is the nth Hermite polynomial H n (x) = (−1) n (exp x 2 ) d n dx n exp(−x 2 ). Consequently φ 0 (x) = π −1/4 exp(−x 2 /2) and 0 = 1/2. The integral kernel of exp(−tĤ 0 ) is a symmetric function of space variables: and stands for a transition density of the diffusion-type process with killing, with a well defined probabilistic interpretation including e.g. the survival time and its exponential decay in the large time asymptotic, see references [30,56]. The kernel induces the semigroup propagation: with a conspicuously time-independent ground state contribution. We note that any suit- can be met only if α 0 = 1. Hence, a proper form for Ψ(x) of interest is In principle, one can work with any Ψ ∈ L 2 (R). The form (A8) can be introduced, if Ψ is not orthogonal to φ 0 . Then, (Ψ, φ 0 ) = α 0 = 0 entails the replacement of Ψ by (1/α 0 )Ψ, which does the job.
We can give a concise summary of of relationships between transition densities k 0 (t, y, x), k(t, y, x) and p(t, y, x), introduced in subsection 2.2. We begin from the killed process, for which the integral kernel of exp(−tĤ 0 ), withĤ 0 = (1/2)(−Δ + x 2 ), reads (the order of space variables is important): Note a conspicuous presence of the time-dependent factor exp(−t/2), comprising the contribution from the lowest eigenvalue 1/2 ofĤ 0 , which is directly responsible for the exponential decay of the transition pdf.
The transition probability density of the relaxation process governed by equations (1), (2), (12), in the present case reads: where φ 0 (x) = ρ 1/2 * (x). Equation (A.10) reproduces the transition density of the Ornstein-Uhlenbeck process in R. Here, k 0 (t, y, x) refers to a process with killing [30,56], while k(t, y, x) refers to the spectrally relaxing process. The presence of the ratio φ 0 (x)/φ 0 (y) is a signature of the (emergent, see reference [30]) Doob-type conditioning which transforms k(t, y, x) into a legitimate (probability measure preserving) transition probability density p(t, y, x) of a Markovian diffusion process with an invariant probability density [30,32,33].

A.4. Interval with absorbing endpoints versus non-absorbing (Dirichlet) infinite well
For the process with killing (absorption) at the boundaries of the interval [−1, 1], the generator of motion is the ordinary Laplacian Δ D with Dirichlet boundary conditions imposed as its domain restriction (L 2 ([−1, 1]) functions, vanishing at x = ±1 [40]). A physically legitimate step is here to analyze the dynamics in terms of a potential function U(x), which mimics the infinite square well enclosure and directly refers to quantum theory routines. Formally, one is inclined to set U(x) = 0 for x ∈ (−1, 1) ⊂ R and U(x) ≡ ∞ for R\[−1, 1] and introduce the operatorĤ 0 = −Δ + U(x) (see [39][40][41]). This may be interpreted as a 'physical' encoding of −Δ D (we recall our introductory mention of the case of a constant potential plus the Dirichlet boundary data imposed upon the operatorĤ, see also [32]).
This well conforms with the general statement [32], that the Schrödinger operator with constant potential and Dirichlet boundary condition corresponds to the Brownian motion that is conditioned never to exit (−1, 1), see also [30].
Recalling the path integration (Feynman-Kac) formula (11), and setting there V = −π 2 /4 (D = 1 case), we realize that − t 0 V(ω(τ ))dτ = +tπ 2 /4 identically. Hence the transition The operatorĤ = −DΔ + x m , with m even (and typically D = 1 or D = 1/2), is known to provide a fairly accurate spectral approximation of the standard infinite well problem with Dirichlet boundaries, whose eigenvalues and eigenfunctions can be reproduced approximately in the form of of 1/m expansions, see e.g. [34,36,37]. Since analytic solutions for moderately sized values of m are unavailable, we have adopted a computer assisted route toward the approximate evaluation of eigenfunctions and eigenvalues ofĤ, with a focus on their convergence properties with the growth of m. Its is based on the Strang splitting method, originally devised in another (Lévy processes and fractional Schrödinger operators) context, see e.g. [11,12].
We reproduce final computer-assisted results, obtained for the case D = 1/2 i.e. for the operatorĤ = − 1 2 Δ + x m . Guided by the 'reconstruction of the dynamics from the ground state' idea, we confine our attention to the m-dependence of the ground state eigenvalue (in the present case it is not equal zero) and the related eigenfunction. The ground state functions are depicted for m = 50, 100, 200, 500, 5000, two latter curves are indistinguishable from the asymptotic cos(πx/2).
The case of x m /m has been covered in a parallel computation, but made comparatively explicit in figure 10, only on the level of eigenvalues. As far as the eigenfunctions are concerned, their qualitative behavior for the potential x m /m is the same as that reported for x m in figure 10.
The presented analysis demonstrates that for large even m, the spectral problem forĤ = −DΔ + x m well approximates the spectral problem for the Dirichlet operator −Δ D on the interval (equivalently, in the infinite well). Since for each m = 2n 2, in the least numerically, we have in hands the ground state function ρ 1/2 * (x), which corresponds to the eigenvalue zero:Ĥρ 1/2 * = (− Δ + V])ρ 1/2 * = 0, we can readily follow the reverse engineering idea (originally devised for Lévy flights) [20,21,42].
Its basic aim has been to reconstruct the Langevin (Fokker-Planck)-driven diffusion process, for which an invariant probability density is given a priori. This is precisely our point of departure, since for a couple of m we have in hands the pertinent ρ * (x).
Our reconstruction procedure goes as follows. Given ρ 1/2 * (x), we take its square. Next we adopt the numerical algorithm allowing to retrieve b(x), through a direct evaluation of b(x) = ∇ ln ρ * = [∇ρ * (x)]/ρ * (x). The results are depicted in figure 12 and show convincingly why (drift spatial behavior) the Langevin-driven Brownian motion in R sets down at an equilibrium pdf. We indeed sequentially approximate the taboo process ( [30]) in the interval, while keeping under control how the asymptotic pdf, drift and Newtonian (Boltzmann, see section 2) potential are approached with the growth of even m. In the interval (alternatively, in the Dirichlet infinite well), the invariant pdf is ρ * (x) = cos 2 (πx/2). The related Fokker-Planck equation is reproduced by evaluating b(x) = ∇ ln ρ * (x) = −∇U(x), with U(x) = −2 ln cos(πx/2) (this corresponds to D = 1/2) playing the role of the Newtonian (and Boltzmann, see section 2) potential ( figure B2).

Appendix C. Schrödinger's boundary data and interpolation problem: arena for conditioned Brownian motions
The previous discussion of diversly confined Brownian motions, has a natural affinity with the classic Schrödinger's boundary data and interpolation problem. Its roots range back to the original paper by Schrödinger [65], the theory being revived by mathematicians and physicists, see e.g. [51][52][53]66] and references therein.
The Schrödinger boundary data problem is known to provide a unique Markovian interpolation between any two strictly positive probability densities, designed to form the input-output statistics data for a random process process bound to run in a finite (observation) time interval.
The key input, if one attempts to reconstruct the pertinent Markovian dynamics, is to select the jointly continuous in space variables, positive and contractive semigroup kernel. Its choice is arbitrary, except for the strict positivity (not a must, but we keep this restriction in the present paper) and continuity demand.
The semigroup dynamics in question we infer from the classic notion of the Schrödinger semigroup exp(−tĤ ), with a proviso that the semigroup generatorĤ actually stands for a legitimate (up to scaled away physical constants) Hamiltonian operator of the form (dimensional constants are scaled away)Ĥ = −Δ + V(x).
Roughly, the essence of the Schrödinger boundary data problem [66], goes as follows. We consider Markovian propagation scenarios running in a finite time interval [0, T ], with the input-output statistics data provided in terms of two strictly positive boundary densities ρ(x, 0) and ρ(x, T ), T > 0 [12,66]. We demand the existence of the bivariate transition probability m (A, B), A, B ⊂ R, admitting the probability density m(x, y) = f (x)k(x, 0, y, T)g(y) ( C 1 ) with marginals R m(x, y)dy = ρ(x, 0) and R m(x, y)dx = ρ(y, T ) that may be constrained to (integrated over) some Borel sets A and B contained in R to yield ρ 0 (A) and ρ(B, T ) as marginals for m(A, B).
In equation (C1), f(x) and g(y) are the a priori unknown strictly positive functions, that need to be deduced from the imposed boundary data (i.e. marginals that are presumed to be known a priori). To this end, in addition to the density boundary data we select a strictly positive, jointly continuous in space variables kernel function k(x, 0, y, T ), which we consider as a restriction (to endpoints of the considered time interval [0, T ]) of a certain strongly continuous dynamical semigroup kernel k(y, s, x, t), 0 s < t T. This assumption will secure the Markov property of the associated interpolating process. Actually, we shall consider time homogeneous processes generated by the semigroup exp[−(t − s)Ĥ ), with the kernel k(t − s, y, x).
Under those circumstances [66] and [12,52,53], once we define functions Here we note the exploitation of the semigroup property in propagation formulas (C2), (C3). Namely, for 0 s < t T we have [66]: θ(x, s) = k(x, s, y, T)g(y) dy = k(x, s, z, t)k(z, t, y, T)g(y) dy dz = k(x, s, z, t)θ(z, t) dz (C7) and likewise for: For a given semigroup exp(−tĤ ) the specific Hamiltonian generatorĤ, the kernel exp(−tĤ )(y, x) = k(t, y, x) and the emerging transition probability density p(t, y, x) of the time homogeneous stochastic process are unique in view of the uniqueness of solutions f(x), g(y) of the Schrödinger boundary data problem [66]. In the case of Markov processes, the knowledge of the transition probability density p(y, s, x, t) (here p(t − s, y, x)) for all intermediate times 0 s < t T suffices for the derivation of all other relevant characteristics of the pertinent random motion.
We have not yet specified any restrictions upon the properties of the potential function V(x), nor its concrete functional form. In the present paper, the potential is expected to be a continuous function and show up definite confining properties, which we set by demanding lim |x|→∞ V(x) = ∞. Then, the Hamiltonian operatorĤ generically admits a positive ground state function ϕ 0 (x) with an isolated bottom eigenvalue λ 0 (typically a fully discrete spectrum is admitted).