Analytical and Numerical Treatments of Conservative Diffusions and the Burgers Equation

The present work is concerned with the study of a generalized Langevin equation and its link to the physical theories of statistical mechanics and scale relativity. It is demonstrated that the form of the coefficients of the Langevin equation depends critically on the assumption of continuity of the reconstructed trajectory. This in turn demands for the fluctuations of the diffusion term to be discontinuous in time. This paper further investigates the connection between the scale-relativistic and stochastic mechanics approaches, respectively, with the study of the Burgers equation, which in this case appears as a stochastic geodesic equation for the drift. By further demanding time reversibility of the drift, the Langevin equation can also describe equivalent quantum-mechanical systems in a path-wise manner. The resulting statistical description obeys the Fokker–Planck equation of the probability density of the differential system, which can be readily estimated from simulations of the random paths. Based on the Fokker–Planck formalism, a new derivation of the transient probability densities is presented. Finally, stochastic simulations are compared to the theoretical results.


Introduction
The Langevin equation was introduced in order to describe the motion of a test particle subjected to a fluctuating force and a viscous drag [1]. Its formulation was later generalized to encompass also other types of systems. The Langevin equation is also fundamental for the stochastic interpretation of Quantum Mechanics (QM) [2] and it also appears, in the form of a geodesic equation, in the scale relativity theory (SR) developed by Nottale [3]. The equation represents a substantial theoretical innovation because it was in fact the first stochastic differential equation. The formal theory of stochastic differential equations was developed much later by the works of Itô and Stratonovich (see, for example, [4] for introduction).
In contrast to the picture of diffusion as an uncorrelated random walk, the theory of dynamical systems makes it possible to treat diffusion as a deterministic dynamical process. There the Langevin dynamics can be also driven by chaotic but deterministic processes [5][6][7]. Emergence of diffusive behavior and Markovian evolution was also addressed by Gillespie [8]. The recent study of Tyran-Kaminska demonstrates that simple diffusion processes can emerge as weak limits of piecewise continuous processes constructed within a totally deterministic framework [7]. This is a finding which lends credence to the widely used techniques of Monte Carlo simulations using pseudo-random number generators. stochastic version of the Newton's law and time reversibility of the process. Interestingly, the form of the stochastic acceleration had to be postulated.
A Lagrangian formulation of stochastic mechanics was achieved by Pavon in complex form [21]. However, the given presentation is far from intuitive. In his treatment, the stochastic Lagrangian is the classical Lagrangian evaluated on a complex-valued velocity field in place of the real-valued classical velocity, while the dynamics is given by a complex-valued stochastic differential equation, similar to the treatment of Nottale. The Lagrangian problem was formulated as a constrained optimization problem, where the dynamics acted as the constraint.

Scale Relativity Theory (SR)
The scale relativity theory extends the principle of relativity also to resolution scales [3,22,23]. The main tenet of the theory is that there is no preferred scale of description of the physical reality. Therefore, a physical phenomenon must be described simultaneously at all admissible scales. While this is consistent with calculus for differentiable signals, the situation changes if non-differentiable models, such as Brownian motion or Mandelbrot's multiplicative cascades [24], are addressed. For these cases, the scale of observation (or resolution) is present irreducibly in the local description of a phenomenon. This led Nottale to postulate the fractality of the underlying mathematical variety (i.e., a pseudo-manifold) describing the observables. It should be noted that, in Nottale's approach, only finite differences are admissible. The scale relativistic approach results in corrections of Hamiltonian mechanics that arise due to the non-differentiability of trajectories, which are treated as virtual paths. Nottale introduces a complex operator that he calls the scale derivative, which acts as a pseudo-derivative (see Section 4 for details).Using this tool, Nottale gives an informal derivation of the Schrödinger equation from the classical Newtonian equation of dynamics, via a quantization procedure that follows from an extension of Einstein's relativity principle called the scale relativity principle.

Stochastic Representation of Trajectories
If one considers the Brownian particle as a subsystem and the surrounding particles as an infinite dimensional thermal reservoir, the Langevin equation precisely models the situation where the subsystem suitably interacts with the thermal reservoir. The type of the effective random force can be identified with a Wiener process, which has continuous but non-differentiable paths almost everywhere. Mathematical descriptions of strongly nonlinear phenomena necessitate the relaxation of the global assumption of differentiability. In contrast, classical physics assumes global smoothness of the signals and continuity of their first two derivatives. Therefore, non-smooth phenomena, such as fractals slip through its conceptual net. This argument can be further elaborated as follows. Consider the measurement of a trajectory in time x(t). Non-differentiability can occur in three scenarios: 1. divergence of the velocity, that is divergence of the difference quotient, 2. oscillatory singularity or 3. difference between forward and backward velocities.
While for scenarios (1) and (2) the velocities (i.e., derivatives) can not be defined mathematically, scenario (3) requires dropping only the assumption of continuity of the resulting velocity. That is, x + (t) = x − (t) at the point of non-differentiability t. A simple example of such behavior is the signal x(t) = |t| around the origin t = 0. While scenario (2) is excluded by the scale relativity theory, scenario (1) leads to scale dependence of the difference quotient. Examples of fractal functions, such as the mathematical Brownian motion paths, are typically of divergent length. This at best can be viewed as a mathematical idealization since in this case the work for moving a particle along its trajectory must be infinite. On the other hand, non-differentiability does not need to occur "everywhere" (i.e., with full Lebesgue measure) on a trajectory. In this case, the trajectory can be almost everywhere differentiable except on a certain dense set of points. Examples of these are the singular functions, such as the Salem-de Rham's functions [25] or the well known Cantor's function. Singular functions have finite lengths, therefore the exerted displacement work is also finite. This makes them promising candidates for conceptualization of non-smooth phenomena in physics.
The relationship between Nelson's and Nottale's approaches can be established in a formal way. For clarity of the argument, we focus on the one-dimensional case. First, let's establish the concept of stochastic embedding of a signal. In the following, we assume that the deterministic signal (i.e., trajectory) will be represented by an equivalence class of stochastic paths having the same expectation as the given deterministic signal. Mathematical notation and preliminaries for the subsequent treatment are presented in Appendix A. A possibly non-differentiable continuous trajectory is represented by a continuous Markov stochastic process evaluated in the virtual space of paths as follows: Definition 1 (Markov Stochastic Embedding). Consider a bounded deterministic signal x(t) on the compact interval T ⊆ R representing time. Define the stochastic embedding S ρ in the probability space (T ⊗ Ω, F , ρ), where ρ is the probability density, as the isomorphism where the random variables sampled at different times t are independent and identically distributed (i. i. d.) and F is a σ-algebra.
Note: the ω-index will be skipped from the notation wherever convenient for clarity. In addition, X t and X(t) will be used interchangeably. Deterministic signals are denoted by the lower case, while the stochastic by upper case letters.
The above definition implicitly assumes that X t ∈ L 1 (T ⊗ Ω, F , ρ) and E ω X(t, ω) < ∞. The name of the embedding is justified by the following Lemma: Lemma 1. The stochastic process under the above definition has the Markov property.
Proof. By construction for fixed t, δ ∈ F The conditional expectation is where ξ ≡ X t+δ is used for notational convenience. However, by independence of the variables ρ(X t+δ , X t ) = ρ(X t+δ )ρ(X t ). Therefore, E ω X t+δ = E ω (X t+δ |X t ). Since δ can be either positive or negative, the claim follows.
Consider the nonlinear problem, where the phase-space trajectory of a system is represented by a Hölder function x(t) (see Appendix A, Definition A2) and t is a real-valued parameter, for example time or curve length. Let us suppose that the continuous temporal evolution of a differential system can be represented by a generalization of the Langevin equation of the form where a(x, t) and B are bounded and measurable functions of the co-ordinates and furthermore a(x, t) is continuous in both x and t. That is, for all , such that 0 ≤ ≤ dt This can be recognized as the Hölder growth condition of order β, since a(x, t) is an O ( ) term. The fractional exponent β is treated as a free parameter with value to be determined later.
The type of admissible functions coupled to the fractional exponent depends critically on the assumption of continuity of the reconstructed trajectory. This in turn demands for the fluctuations of the fractional term to be discontinuous. The proof technique is introduced in [26], while the argument is similar to the one presented by Gillespie [8].
Without loss of generality, set a = 0. Let x t+ = x t + B(x t , t) + O β and |∆ x| ≤ K β . Fix the interval [t, t + ] and choose a partition of points P = {t k = t + k/N} Therefore, by induction If we suppose that B is continuous in x, implying also continuity in t, after taking supremum limit on both sides lim sup Therefore, either β = 1 (which is forbidden by hypothesis) or else B = 0 so that B(x, t) must oscillate from point to point if β < 1. Then, let's denote the set χ β := {B(x t , t) = 0}.
The argument demonstrates that so-defined set is totally disconnected in the topology of the real line [26]. This allows for the choice of the algebra F , since we can demand that Ω ⊆ χ β has for elements the semi-open intervals [τ i , τ j ), τ i,j ∈ χ β . Furthermore, the initial system in Equation (1) is equivalent to the finite existence of the fractional velocity B(x, t) = υ β + x(τ i ) = 0, since the differential system can be recognized as fractional Taylor series [26]. In other words, the events in the probability space are the observations of non-vanishing values of the fractional velocity of the signal.
From now on, let P τ ≡ P ⊆ F . Without loss of generality, suppose that O N 1−β β ≤ 1. The stochastic representation x t → (X t (ω), ρ) is such that Therefore, we demand that B(X t , t) is F -measurable and L 2 (Ω, T) as a technical condition. By the Hölder condition, |x t k − x t k−1 | ≤ K k β for some set of constants K k . Then, by transfer, Therefore, E B(X t k , t k ) exists and is bounded. By the same argument, Then, we proceed by induction. Let K s = sup i K 2 i from the above partition: Therefore, for the embedded variable, Since ∆ i X∆ j X are independent by Lemma 1, E∆ i X∆ j X = E∆ i XE∆ j X ≤ K s . Therefore, The argument can be specialized to β = 1/2 where Var[∆ X] ≤ K s 2β . Therefore, the variance exists ∀N and the Central Limit Theorem holds. Since by Lemma 1 the process is Markovian, it must follow that in limit N → ∞ the random process is Wiener. Now suppose that a = 0. Then, since a(x) is continuous of bounded variation (BVC, see Appendix A), then a.e., with σ 2 existing by the previous argument. Therefore, 2β N 1−2β by the same argument as in the previous case. Therefore, for β = 1/2, the limit of the random process is Wiener.
Let us denote the limit Wiener process by W t . Using the stationarity and self-similarity of the increments ∆ + W t = √ N(0, 1), where N(0, 1) is a standard Gaussian random variable. Therefore, for β = 1/2, the velocity can be regularized to a finite value if we take the expectation. That is, The estimate holds a.s. since P(W t = 0) = 0, where P denotes probability. Finally, there is a function b(X, t), such that b(X, t)ξ = B(X, t), ξ ∼ N(0, 1). This follows directly from the axiom of choice, since we can always choose ξ = 1. Therefore, the last equation can be treated as a definition of b(X, t).
In summary, the following theorem can be formulated: is continuous in both x and t and B(x, t) ∈ L 2 (Ω, T) is bounded but discontinuous. Furthermore, let χ β be the set of change (Definition A6) of f [T].
Such embedding can be also called a consistent stochastic embedding. This theorem allows for Nelson's characterization of the Langevin diffusion process.

Nelson's Characterization
The Langevin equation can describe equivalent quantum-mechanical systems in a path-wise manner. These are the so-called conservative diffusions of Carlen [27]. The existence of so-conceived QM particle paths was proven under certain reasonable conditions [27]. Starting from the generalized Langevin equation, the argument can be specialized to a Wiener driving process, which can be handled using the apparatus of Itô calculus. Consider the stochastic differential equation with continuous drift and diffusion coefficients where a(X, t) and b(X, t) are smooth functions of the co-ordinates and dW t are the increments of a Wiener process dW dt ∼ N(0, dt) adapted to the past filtration F t>0 -i.e., starting from the initial state. Let EX t = x(t). Following Nelson [2], the forward and backward and drift, respectively diffusion coefficients, can be identified as the averaged velocities [28]: The evolution of the density of the process can be computed from the forward Fokker-Planck equation which can be recognized as a conservation law for the probability current j: Under the finite energy technical condition, there is a backwards process with the same transition density which is adapted to the future filtration F t<T -i.e., starting from the final state. This leads to the anticipative (i.e., anti-Itô or anticipative) stochastic integrals. This process has Fokker-Planck equation Then, it follows that the Nelson's osmotic velocity can be defined from where φ(t) is an arbitrary C 1 function of time as u := 1 2 (a −â) and the current velocity as v := 1 2 (a +â) so that a continuity equation holds for the density Furthermore, Pavon [21] has established that the entropy production over the whole space is Thus, for a Markov diffusion process, for a constant b.

The Complex Velocity Operator in SR and SM Theories
Scale relativity treats velocity only as a difference quotient. This is a necessity due to the assumed non-differentiability of the trajectories. Non-differentiability leads to introduction of two velocity fields-forward and backward, depending on the direction of differentiation in time. These fields are assumed to be finite for small values of the time step dt but they diverge to infinity in the limit dt → 0 in a standard analysis setting. Therefore, such velocity fields can be defined only up to a finite resolution underlying the physical phenomenon under study. The velocity fields are assumed to admit representation of the form of a sum of a "classical part" plus a correction of a resolution-dependent and diverging fractal part. The classical part corresponds to the absolutely continuous part of the trajectory, while the fractal part corresponds to the singular and possibly oscillatory parts. Since, at the level of physical description, there is no way to favor the forward rather than the backward velocity, the description should incorporate them on equal grounds, i.e., forming a bivariate vector field R ⊗ R This bivariate vector field is represented by a complex-valued vector field [29] as v = V − iU ∈ R 3 with components given by U : where V is interpreted as the "classical" velocity and U is a new quasi-velocity quantity (i.e., the osmotic velocity in the terminology of Nelson). Under these assumptions, Nottale introduces a complexified material derivative, which is a pseudo-differential operator acting on scalar functions as where σ is a constant, quantifying the effect of changing the resolution scale.
Stochastic mechanics allows for a similar treatment of the complex of forward and backward diffusions. The drift, resp. diffusion coefficients can be further embedded in a complex space as proposed by Pavon [21]: so that the diffusion process becomes complex. It follows that In the case when b =b, Therefore, we can designate a new complex stochastic variable Because of its double adaptation, Z t retains its local martingale properties: that is, EZ t = 0. In this case, notably Var Z t = 0, but E|Z t | 2 = 1, so that finally, Therefore, a formal Itô differential can be introduced in exactly the same way with quadratic variation dX 2 = −ib 2 dt. Therefore, in components, which generalize to in three dimensions [28]. It is apparent that both theories share an identical algebraical structure, while SM can be considered as a stochastic representation of SR.

Remark 1.
Conceptually, the forward process can be interpreted as a prediction, while the backward process can be interpreted as a retrodiction. Note that, in the complex formulation of Pavon, the real part of the driving process Z t corresponds to the forward (i.e., adapted to the past) process, while the imaginary part corresponds to the backward (i.e., adapted to the future) process. This is of course one of infinitely many choices, since the complex factor in the diffusion coefficient is a root of unity and hence represents a rotation in the complex plane.
The martingale property of the complex Wiener process conceptually means that the knowledge of the past and future of the process do not bias the outcome at the present time (i.e., at measurement). Note that the mapping is invertible since From these formulas, it is apparent that the real part, or respectively the imaginary part of the resulting process do not have separate meanings, as they mix the predictive process with the retrodictive process. To illustrate the point, suppose that F = F r + iF i and a = a r + ia i and the original process dX is transformed as F(X ). Then, a straightforward calculation gives

The Real Stochastic Geodesic Equations
The appearance of the Wiener process entails the application of the fundamental Itô Lemma for the forward (i.e., adapted to the past, plus sign) or the backward processes (i.e., adapted to the future, minus sign), respectively. In differential notation, it reads where dX 2 = b 2 dt is the quadratic variation of the process. It can be seen that in this case the (forward) differential operator d acts as a material derivative.
The term geodesic will be interpreted as a solution of a variational problem [30,31]. A brief treatment is given in Appendix B. By application of Itô's Lemma, the forward geodesic equation can be obtained as: This can be recognized as a Burgers equation with negative kinematic viscosity for the drift field [13].
The backward geodesic equation follows from the application of the Itô's lemma for the anticipative process This can be recognized as a Burgers equation with positive kinematic viscosity for the drift field. The solution of the Burgers equation is well known and can be given by the convolution integrals (Equation (44)) for the case of positive viscosity [13]. The case about the negative viscosity can not be easily solved using Fourier transform. Therefore, a different solution technique will be pursued. Time-varying solutions will be constructed from topological deformations of the stationary solutions.
In QM applications, b =h/2m. Normalization b = 1 will be assumed further in most cases to simplify calculations.

Path-Wise Separable Solutions
In the first instance, one can solve the geodesic equation by supposing separability. By making the ansatz a(x, t) = f (x)g(t), we arrive at the equation: This has the unique solution The resulting Itô equation can be formulated as The stochastic differential equation for the drift is therefore which can be integrated exactly in Itô's sense as Therefore, where T is the stopping time. Therefore, an exact numerical quadrature can be performed (Figure 1   The corresponding density can be obtained from the Fokker-Planck equation .
It should be noted that, under time reversion, we arrive at the same solution, which however leads to a different Fokker-Planck equation which can be recognized as a Brownian bridge. The entropy of this density can be calculated as

Stationary Drift Fields
For time-homogeneous diffusion, the geodesic equation can be brought into the form which can be integrated once to give The integration constant E can be identified with the energy. The resulting first order ordinary differential equation (ODE) can be solved as The solution for E > 0 was identified by Herman [32]. By translation, invariance of the coordinates, c = 0 is admissible. This observation will be used further for the transient solution. The link between the two solutions can be established as follows. Note that is also a solution. Then, which is the second solution.
The expectation of the trajectory can be obtained by solving the ODE In accordance with so-developed theory for c = 0, so that in the same way The backward geodesic equation by the same method leads to

Stationary Density Solutions
The stationary density ρ(y) is a solution of the Fokker-Planck (i.e., forward Kolmogorov) equations parametrized by E: The case E > 0 leads to The case E = 0 leads to with a stationary solution ρ = |x|, which can be valid on a bounded domain.

Transient Drift Fields
The solution of the Burgers equation is well known and can be given by the convolution integrals (Equation (44)) for the case of positive viscosity [13]. The case of negative viscosity emerging here is more challenging and it will be solved by a deformation of the stationary solution, so that in limit the stationary solution is recovered: lim The solution is sought in the form (neglecting scale factors) which results in a linear ODE for the unknown function f (t): By variation of the parameters, the solution for a(t, x) is given as where the constant E represents an energy scale and k is an arbitrary constant. We can assume normalization, for example k = ± √ E, such that a(t, π 2E ) = ±1. Plots are presented in Figure 2. The transformed Itô drift equation for k = 1 reads It can be further noticed that rescaling in a pair of new variables x = √ Ex, t = Et leaves the ratio invariant so that z becomes a similarity variable. Furthermore, a formal forward Kolmogorov equation can be written in the y = a(t, x) variable with E = 1 as however its solution is challenging due to its mixed nonlinearity and will not be attempted here. Nevertheless, the analysis presented so far assures that asymptotically ρ can be obtained as a solution of the stationary equation.
The backward geodesic equation leads to the following solution :

Asymptotic Density Solutions
The forward drift itself a ≡ y (symbol changed) obeys the transformed stochastic differential equations The density ρ is a solution of the forward Kolmogorov equations parametrized by E: The solutions can be obtained using the Laplace transform L s f (t) →f (s). In this way, the partial differential equation can be transformed into an ODE for the Laplace variable: To obtain the Green's function, we take homogeneous initial conditions a.e. The solutions in the time domain can be obtained by the inverse Laplace transformation: ρ(s, y) = e − √ 2s/E−1 arctan(y) (y 2 + 1) In the position space, the solution can be obtained using Grisanov's theorem [4]: The second equation is not acceptable from a physical point of view since lim t→∞ ρ(t, x) diverges. In the same way for the backward drift, with solutions in the drift space and in position space respectively. This is acceptable from a physical point of view since lim t→∞ ρ(t, x) = 0, which is a correct asymptotic behavior.

The Complex Stochastic Geodesic Equations
The complexification removes the restriction of positive definiteness of the E parameter so that the substitution t → ±Et becomes admissible by an appropriate cut along the complex plane.
In a similar way, for the complex case, we have which, under substitution, y = tanh x leads to By the same methods as used above, the asymptotic density for the drift variable can be obtained as For the resulting density in the position space, it can be calculated that In a similar way, for the other solution, which under substitution y = tan √ Ex leads to the drift equation The drift density can be readily obtained as In the position space, the density is of the form In either case, the densities asymptotically approach zero.

Real-Valued and Complex Cole-Hopf Transformations
The Burgers equation can be linearized by the Cole-Hopf transformation [33,34]. This mapping transforms the nonlinear Burgers equation into the linear heat conduction equation in the following way. Let Substitution into Equation (11) leads to This can be recognized as ∂ ∂x ∂ 2 ∂x 2 u = 0, which is equivalent to a solution of the equation It should be noted that if instead of the forward development (i.e., prediction) one takes the backward development (i.e., retrodiction), the usual form of the Burgers equation is recovered. This corresponds to the anticipative Wiener process, which is subject to the anticipative Itô calculus [17,35]: In this case, the usual general solution can be revealed a(x, t) = ∂ ∂x where ν = 1/2 is the viscosity coefficient.
In the complex case, starting from the generalized Itô differential, the complex velocity field becomes The geodesic equation reads Therefore, by the martingale property, this is equivalent to which can be recognized as a generalized Burgers equation with imaginary kinematic viscosity coefficient. Applying the complex Cole-Hopf transformation as [36] V = −i ∂ ∂x log U, −π < arg U < π and specializing to b = 1 leads to which can be recognized as a gradient The last equation is equivalent to the solution of the free Schrödinger equation. On the other hand, the diffusion part is simply This corresponds with the arguments given in [37] that the coefficient of the stochastic noise should be purely imaginary. Calculations can be reproduced in the computer algebra system Maxima [38].

Numerical Results
The different types of solutions of the stochastic geodesic equation were simulated using the Euler-Maruyama algorithm. Simulations were performed in Matlab. An example of a simulation script is given in Appendix C.

Exact Simulations
The exact simulations of the separable process were compared to simulations computed by the Euler-Maruyama algorithm. Achieved correlation was 1.0 while the mean squared error was on the order of 1 × 10 −8 .The comparison is presented in Figure 1. The empirical transition density was computed from N s = 1000 simulations and correlated to the graph of Equation (16). The theoretical density was computed from Equation (16) for stopping time T = 10. Achieved correlation was 0.9976.

Free Diffusion
The normalized asymptotic transient density of the free particle distribution can be recognized as the Rayleigh's distribution ( Figure 3)

Particle in a Box
The third simulated case comprised a freely diffusing particle in a square potential well of size 2L. The approach was based on Hermann [32]. Individual trajectories were simulated according to the fundamental equation using the scheme of Euler-Maruyama: where ∆W n ∼ N(0, 1). Restarting boundary conditions were used for the simulations to avoid distortions of the distribution. That is, if a simulated particle crossed the boundaries its position was reset to its original position.
The initial particle positions were sampled from a uniform distribution between −L and L. The theoretical density for the particle in a box case is given by Results are based on N = 10000 points in N s = 1000 simulations. The empirical pdf is estimated from n = log(N s N) 2 bins. Pearson's correlations are given as insets: B -r = 0.9939, D -r = 0.9871. For both cases, the numerical precision correlates excellently with the analytical solutions ( Figure 4).

Discussion
This work was motivated in part by the premise that inherently nonlinear phenomena need development of novel mathematical tools for their description. The relaxation of the differentiability assumption opens new avenues in describing physical phenomena, as demonstrated by SM and SR, but also challenges existing mathematical methods, which are developed for smooth signals [2,3]. While this description can be achieved also by fractional differ-integrals, or by multi-scale approaches [39], the present work focused on a local description. The reason for this choice is that locality provides a direct way of physical interpretation of the obtained results. In this regard, Hölderian functions can be used as building blocks of such strongly nonlinear models, which give rise to singular [24,40] or non-differentiable models.
The second motivation of the present work was to investigate the potential of stochastic methods for simulations of quantum-mechanical and convection-diffusive systems. While the usual presentation of the stochastic mechanics typically used the Schrödinger equation as a solution device and paths were constructed from solutions of the Schrödinger equation, this is not necessary. McClendon and Rabitz simulated several quantum systems using the differential equations of Nelson's stochastic quantization as a starting point [41]. In the framework of scale relativity, Herman [32] and later Al Rashid et al. [42] simulated QM particle in a box using the Langevin equations. Later, Al-Rashid et al. [43] simulated the quantum harmonic oscillator extending Herman's approach. The approach presented here can be used as an alternative to numerical solutions of the Schrödinger equation. In this scenario, the density of the solution can be sampled from Monte Carlo simulations as demonstrated. Presented numerical approaches can be used, for example, for simulations of nanoparticles or quantum dots, which are mesoscopic objects and are expected to have properties intermediate between macroscopic and quantum systems [44]. This can be of interest, for example in sedimentation studies, where Langevin dynamics was proposed [45]. In principle, presented results can be extended towards asynchronous simulations using the Gillespie's algorithm [8]. This can be achieved using time steps distributed exponentially.
Obtained results can be also discussed in view of the fluctuation-dissipation relationships. The fluctuation-dissipation theorem relates the linear response relaxation of a system from a non-equilibrium state to the properties of fluctuations in equilibrium. This is an exact result in the case of the Ornstein-Uhlembeck process, where the drift term is linear. The geodesic treatment in the present work provides a different relationship between the drift and a(x, t) and diffusivity b(x, t). In the small perturbation regime around the equilibrium the geodesic process x eq (t) can be approximated by an Ornstein-Uhlembeck process for the fluctuation term (ξ = δx(t)), therefore an appropriate fluctuation-dissipation theorem can be formulated assuming that equipartition also holds.
A fact that is not fully addressed by both stochastic mechanics and scale relativity is why do the theories work only for (box) fractal dimension 2 of the paths. While Nottale gives an heuristic argument and claims that the prescription of a Wiener process may be generalized, he does not proceed to rigorously develop the argument. On the other hand, the stochastic mechanics fixes from the start the Wiener process as a driving noise. While this may look plausible in view of the traditions in the treatment of Brownian motion, it is a choice that should be justified as nowadays anomalous types of diffusion dynamics are also recognized and systematically investigated (overview in [46]). The answer to this question can be given more easily by an approach inspired by Nottale and is partially given by the argument given by Gillepsie [8]. The original argument in [8] contains an explicit assumption of existence of the second moment of the distribution, which amounts to assuming Hölder continuity of order 1/2 as demonstrated here in Theorem 1. The theorem also corresponds to the result established for fractal interpolation computed via a chaos game where the limit random distribution has been identified with the Gaussian distribution [47,48].
While in SR particle 'trajectories' are considered to be only virtual, SM and the original formulation of Bohm's quantum mechanics treat them as physically real. It is noteworthy that recently Flack and Hiley [49] demonstrated that that a Bohm 'trajectory' is the average of an ensemble of actual individual stochastic Feynman paths. This is in line with the treatment of the problem by the stochastic mechanics and scale relativity and promotes the view that Bohm's quantum mechanics is a mean field theory of the stochastic mechanics.
Funding: The work has been supported in part by grants from Research Fund-Flanders (FWO), contract number VS.097.16N, Fetzer Franklin Memorial Trust and the NanoStreeM H2020 project, contract number 688194.

Definition A2.
We say that f is of (point-wise) Hölder class H β if for a given x there exist two positive constants C, δ ∈ R that for an arbitrary y ∈ Dom[ f ] and given |x − y| ≤ δ fulfill the inequality | f (x) − f (y)| ≤ C|x − y| β , where | · | denotes the norm of the argument. Definition A3. Define the parametrized difference operators acting on a function f (x) as The first one we refer to as forward difference operator, the second one we refer to as backward difference operator.
Definition A4. Define Fractional Variation operators of order 0 ≤ β ≤ 1 as This section follows the presentation given recently in [50].
Definition A5 (Fractional order velocity). Define the fractional velocity of fractional order β as the limit where 0 < β ≤ 1 are real parameters and f (x) is real-valued function. A function for which at least one of υ β ± f (x) exists finitely will be called β-differentiable at the point x.
In the above definition, we do not require upfront equality of left and right β-velocities. This amounts to not demanding continuity of the β-velocities in advance. Instead, continuity is a property, which is fulfilled under certain conditions. Definition A6. The set of points where the fractional velocity exists finitely and υ β ± f (x) = 0 will be denoted as the set of change χ β ± ( f ) := x : υ β ± f (x) = 0 .
Since the set of change χ α + ( f ) is totally disconnected [26], some of the useful properties of ordinary derivatives, notably the continuity and the semi-group composition property, are lost. Definition A7. β-Regularized derivative of a function is defined as: We will require as usual that the forward and backward regularized derivatives be equal for a uniformly continuous function.
In this section, we assume that the functions are BVC in the neighborhood of the point of interest. Under this assumption, we have