# Short-time expansion of one-dimensional Fokker-Planck equations with heterogeneous diffusion

Tom Dupont, Stefano Giordano, Fabrizio Cleri, Ralf Blossey

► To cite this version:

Tom Dupont, Stefano Giordano, Fabrizio Cleri, Ralf Blossey. Short-time expansion of one-dimensional Fokker-Planck equations with heterogeneous diffusion. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2024, 109 (6), pp.064106. 10.1103/PhysRevE.109.064106. [https://arxiv.org/abs/2006.06410] (https://arxiv.org/abs/2006.06410)

HAL Id: hal-04614378

https://hal.science/hal-04614378

Submitted on 17 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Short-time expansion of one-dimensional Fokker-Planck equations with heterogeneous diffusion

Tom Dupont\* and Stefano Giordano†

University of Lille, CNRS, Centrale Lille, Univ. Polytechnique Hauts-de-France,
UMR 8520 - IEMN - Institut d'Électronique, de Microélectronique et de Nanotechnologie, F-59000 Lille, France

Fabrizio Cleri‡

University of Lille, Institut d'Électronique, de Microélectronique et de Nanotechnologie
(IEMN CNRS UMR8520) and Departement de Physique, F-59652 Villeneuve d'Ascq, France

Ralf Blossey§

University of Lille, Unité de Glycobiologie Structurale et Fonctionnelle (UGSF), CNRS UMR8576, F-59000 Lille, France

(Dated: 29 mai 2024)

We formulate a short-time expansion for one-dimensional Fokker-Planck equations with spatially dependent diffusion coefficients, derived from stochastic processes with Gaussian white noise, for general values of the discretization parameter $0 \le \alpha \le 1$ of the stochastic integral. The kernel of the Fokker-Planck equation (the propagator) can be expressed as a product of a singular and a regular term. While the singular term can be given in closed form, the regular term can be computed from a Taylor expansion whose coefficients obey simple ordinary differential equations. We illustrate the application of our approach with examples taken from statistical physics and biophysics. Further, we show how our formalism allows to define a class of stochastic equations which can be treated exactly. The convergence of the expansion cannot be guaranteed independently from the discretization parameter $\alpha$.

## I. INTRODUCTION

The Fokker-Planck (FP) equation is one of the most important computational tools in statistical mechanics to obtain an approximation of the probability distribution in time and space of a generic system including a stochastic driving force. As such, the FP equation has found numerous applications in widely disparate domains, ranging from plasma physics, biophysics, physical chemistry, finance, and others [1–3]. Interestingly, the underlying theory of this linear parabolic partial differential equation is still not yet fully developed, and continues to hold surprises of physical relevance. Recent developments e.g. concern spatially varying diffusivity adopted to describe anomalous diffusion [4–12]. Among the many statistical properties of heterogeneous diffusion processes studied are in particular non-stationarity and non-ergodicity; see e.g. results on the non-existence of normalizable stationary solutions, or equilibrium probability distributions, which can be tackled by the methods of infinite ergodicity theory [13–16]. In particular these latter results show a sensitive dependence on the discretization parameter $\alpha$ of the stochastic integrals in the Langevin equation underlying the FP equation. The two most commonly chosen discretizations are the Itô and Fisk-Stratonovich, using values of $\alpha=0$ or $1/2$, respectively.

These recent results raise the question whether a similar dependence on the discretization parameter $\alpha$ exists when considering the short-time behaviour of the FP equation, as it can be deduced from systematic expansions in this limit. Short-time expansions have been considered before in the literature

[17–25], however, to the best of our knowledge, not with this specific question in mind.

Apart from this important formal aspect, our work was also originally motivated by a specific biophysical problem, namely statistical physics models for molecular motors, in which FP equations with heterogeneous diffusion terms arise [26, 27]. In this context we are particularly interested in chromatin remodelers, molecular motors and enzymes that move and remove nucleosomes in eukaryotic DNA [28]. We will treat this problem as one of our applications of our approach in this paper.

It is interesting to point out that multiplicative noise has played an essential role in modeling a number of physical and biological phenomena among which we can mention the transmission of signals in neuron models [29, 30], the phenotypic variability and gene expression [31, 32], the stochastic thermodynamics of holonomic systems [33–35], the ballistic-to-diffusive transition of the heat propagation [36, 37], and the statistical theory of turbulence [38, 39]. Given the widespread use of models based on multiplicative noise and heterogeneous diffusion, it is important to have exact or approximate techniques to solve the corresponding Langevin and Fokker-Planck equations and thus trace the statistical properties of the stochastic processes involved. The approach proposed in this paper, concerning a development of the propagator for short timescales, thus fits into this context and is complementary to all other solution techniques proposed in the literature.

Our paper is organized as follows. We group our material into three sections. Section II begins with a precise formulation of the mathematical problem we discuss and the subsequent derivation of the equations of our formal method. Section III presents the explicit discussion of four exemplary stochastic processes, which illustrates the application of our method with concrete examples. In Section IV we turn our discussion around by exploiting the properties of our expan-

<footnote>\* tom.dupont@centrale.centralelille.fr</footnote>

<footnote>† stefano.giordano@univ-lille.fr</footnote>

<footnote>‡ fabrizio.cleri@univ-lille.fr</footnote>

<footnote>§ ralf.blossey@univ-lille.fr</footnote>
# Short-time expansion of one-dimensional Fokker-Planck equations with heterogeneous diffusion

Tom Dupont\* and Stefano Giordano†

University of Lille, CNRS, Centrale Lille, Univ. Polytechnique Hauts-de-France,
UMR 8520 - IEMN - Institut d'Électronique, de Microélectronique et de Nanotechnologie, F-59000 Lille, France

Fabrizio Cleri‡

University of Lille, Institut d'Électronique, de Microélectronique et de Nanotechnologie
(IEMN CNRS UMR8520) and Departement de Physique, F-59652 Villeneuve d'Ascq, France

Ralf Blossey§

University of Lille, Unité de Glycobiologie Structurale et Fonctionnelle (UGSF), CNRS UMR8576, F-59000 Lille, France

(Dated: 29 mai 2024)

We formulate a short-time expansion for one-dimensional Fokker-Planck equations with spatially dependent diffusion coefficients, derived from stochastic processes with Gaussian white noise, for general values of the discretization parameter $0 \le \alpha \le 1$ of the stochastic integral. The kernel of the Fokker-Planck equation (the propagator) can be expressed as a product of a singular and a regular term. While the singular term can be given in closed form, the regular term can be computed from a Taylor expansion whose coefficients obey simple ordinary differential equations. We illustrate the application of our approach with examples taken from statistical physics and biophysics. Further, we show how our formalism allows to define a class of stochastic equations which can be treated exactly. The convergence of the expansion cannot be guaranteed independently from the discretization parameter $\alpha$.

## I. INTRODUCTION

The Fokker-Planck (FP) equation is one of the most important computational tools in statistical mechanics to obtain an approximation of the probability distribution in time and space of a generic system including a stochastic driving force. As such, the FP equation has found numerous applications in widely disparate domains, ranging from plasma physics, biophysics, physical chemistry, finance, and others [1–3]. Interestingly, the underlying theory of this linear parabolic partial differential equation is still not yet fully developed, and continues to hold surprises of physical relevance. Recent developments e.g. concern spatially varying diffusivity adopted to describe anomalous diffusion [4–12]. Among the many statistical properties of heterogeneous diffusion processes studied are in particular non-stationarity and non-ergodicity; see e.g. results on the non-existence of normalizable stationary solutions, or equilibrium probability distributions, which can be tackled by the methods of infinite ergodicity theory [13–16]. In particular these latter results show a sensitive dependence on the discretization parameter $\alpha$ of the stochastic integrals in the Langevin equation underlying the FP equation. The two most commonly chosen discretizations are the Itô and Fisk-Stratonovich, using values of $\alpha=0$ or $1/2$, respectively.

These recent results raise the question whether a similar dependence on the discretization parameter $\alpha$ exists when considering the short-time behaviour of the FP equation, as it can be deduced from systematic expansions in this limit. Short-time expansions have been considered before in the literature

[17–25], however, to the best of our knowledge, not with this specific question in mind.

Apart from this important formal aspect, our work was also originally motivated by a specific biophysical problem, namely statistical physics models for molecular motors, in which FP equations with heterogeneous diffusion terms arise [26, 27]. In this context we are particularly interested in chromatin remodelers, molecular motors and enzymes that move and remove nucleosomes in eukaryotic DNA [28]. We will treat this problem as one of our applications of our approach in this paper.

It is interesting to point out that multiplicative noise has played an essential role in modeling a number of physical and biological phenomena among which we can mention the transmission of signals in neuron models [29, 30], the phenotypic variability and gene expression [31, 32], the stochastic thermodynamics of holonomic systems [33–35], the ballistic-to-diffusive transition of the heat propagation [36, 37], and the statistical theory of turbulence [38, 39]. Given the widespread use of models based on multiplicative noise and heterogeneous diffusion, it is important to have exact or approximate techniques to solve the corresponding Langevin and Fokker-Planck equations and thus trace the statistical properties of the stochastic processes involved. The approach proposed in this paper, concerning a development of the propagator for short timescales, thus fits into this context and is complementary to all other solution techniques proposed in the literature.

Our paper is organized as follows. We group our material into three sections. Section II begins with a precise formulation of the mathematical problem we discuss and the subsequent derivation of the equations of our formal method. Section III presents the explicit discussion of four exemplary stochastic processes, which illustrates the application of our method with concrete examples. In Section IV we turn our discussion around by exploiting the properties of our expan-

<footnote>\* tom.dupont@centrale.centralelille.fr</footnote>

<footnote>† stefano.giordano@univ-lille.fr</footnote>

<footnote>‡ fabrizio.cleri@univ-lille.fr</footnote>

<footnote>§ ralf.blossey@univ-lille.fr</footnote>
sion in order to find stochastic processes that are explicitly
solvable. We conclude in Section V with an outlook on fur-
ther possible developments.

## II. THE SHORT-RANGE EXPANSION OF THE FOKKER-PLANCK EQUATION WITH HETEROGENEOUS DIFFUSION

In this Section, will develop the formalism to study the short-time behavior of the one-dimensional Fokker-Planck equation. We initially state the problem under investigation and then we subdivide the development into two steps : (i) firstly, we obtain the singular short-time part of the propagator, and (ii) secondly we get the differential equations for the regular larger-time part. Finally, we obtain the formal solution of the hierarchy of differential equations describing the Taylor coefficients of the regular part of the propagator, particularly useful for the applications.

### A. Mathematical formulation of the problem

The starting point of our study is the Langevin equation, a stochastic differential equation of the form

$$ \frac{dx}{dt} = h(x,t) + g(x,t)\xi(t), \quad (1) $$

where $x(t)$ is the stochastic process under investigation, and $h$ and $g$ are regular functions. Without limiting generality, we can always assume that $g(x,t) \ge 0$ for any $x \in \mathbb{R}$ and $t \ge t_0$.

Moreover, $\xi(t)$ in Eq. (1) is a Gaussian white noise with average value $\mathbb{E}(\xi(t)) = 0$, and correlation $\mathbb{E}(\xi(t)\xi(t')) = 2\delta(t-t')$, where $\delta(t)$ is the Dirac delta function. Throughout this work, we always consider processes defined on the entire real axis, without any reflecting or absorbing boundary conditions.

For the process $x(t)$, we can introduce the probability density $W(x,t)$. It means that for any real numbers $x$ and $y$ (such that $y \le x$) we determine the probability that $y \le x(t) \le x$ as $\mathbb{P}(y \le x(t) \le x) = \int_y^x W(z,t)dz$, or equivalently, $W(x,t) = \frac{\partial \mathbb{P}(x(t) \le x)}{\partial x}$.

In order to completely define the meaning of Eq. (1), we must introduce a parameter $\alpha \in [0,1]$, identifying the integral stochastic interpretation adopted [1–3, 40, 41]. Indeed, the value of $\alpha$ defines the position of the point at which we calculate any integrated function in the small intervals of the adopted Riemann sum. This integration rule includes $\alpha = 0$ (Itô integral) [42], $\alpha = \frac{1}{2}$ (Fisk-Stratonovich integral) [43, 44], and $\alpha = 1$ (Hänggi-Klimontovich integral) [45, 46] (see also Ref. [47] for some comparison of these different choices).

The probability density $W(x,t)$ is the solution of the Fokker-Planck equation [1–3, 48–50]

$$ \frac{\partial W}{\partial t} = -\frac{\partial}{\partial x} \left[ \left( h + 2\alpha g \frac{\partial g}{\partial x} \right) W \right] + \frac{\partial^2}{\partial x^2} (g^2 W), \quad (2) $$

which can be also written in the equivalent form

$$ \frac{\partial W}{\partial t} = \frac{\partial}{\partial x} \left[ -hW + g^{2\alpha} \frac{\partial}{\partial x} \left( g^{2(1-\alpha)} W \right) \right]. \quad (3) $$

The right hand side of this Fokker-Planck equation can be expressed more explicitly as linear combination of $W$, $\frac{\partial W}{\partial x}$, and $\frac{\partial^2 W}{\partial x^2}$. We obtain

$$ \frac{\partial W}{\partial t} = -fW - \ell \frac{\partial W}{\partial x} + g^2 \frac{\partial^2 W}{\partial x^2} \quad (4) $$

where we introduced the regular functions

$$ f = \frac{\partial h}{\partial x} + 2(\alpha - 1) \left( \frac{\partial g}{\partial x} \right)^2 + 2(\alpha - 1) g \frac{\partial^2 g}{\partial x^2}, \quad (5) $$

$$ \ell = h + 2(\alpha - 2)g \frac{\partial g}{\partial x}. \quad (6) $$

These expressions show that the Fokker-Planck equation, and thus the density $W$, are not affected by the sign of the function $g$, which we assumed as always positive.

To solve the Fokker-Planck equation, we have to specify the initial condition $W(x,t_0) = W_0(x) \forall x \in \mathbb{R}$, at a fixed initial time $t_0$. In particular, a special solution is given by the kernel or propagator, defined by the initial condition $W(x,t_0) = \delta(x - y)$, where $y$ is an arbitrary initial value. It means that the initial value is deterministically known. This solution is referred to as $W(x,t) = K(x,t;y,t_0)$ and it can be obtained by the partial differential problem

$$ \frac{\partial K}{\partial t} = -fK - \ell \frac{\partial K}{\partial x} + g^2 \frac{\partial^2 K}{\partial x^2} \quad (7) $$

$$ K(x, t_0; y, t_0) = \delta(x - y). \quad (8) $$

With a little abuse of language, we can say that $K$ is the probability density of observing the value $x$ at time $t$, given the value $y$ at time $t_0$.

The propagator is important since any initial condition, such as $W(x,t_0) = W_0(x) \forall x \in \mathbb{R}$, can be handled by the relation

$$ W(x,t) = \int_{-\infty}^{+\infty} K(x,t;y,t_0)W_0(y)dy, \quad (9) $$

as a consequence of the linearity of the Fokker-Planck equation.

The objective of this work is to find a good approximation of the propagator $K$ for short times, and this will be organized as follows. The short-time approximation obviously contains the singularity induced by the initial condition represented by the delta function $\delta(x - y)$. It means that we look for $K_0$ such that $K(x,t;y,t_0) \sim K_0(x,t;y,t_0)$. Successively, for larger values of time, we can introduce a correction function $F$ such that $K(x,t;y,t_0) = K_0(x,t;y,t_0)F(x,t;y,t_0)$. Since the effects of the initial singularity are all contained in $K_0$, it follows that $F$ is regular. We will obtain a partial differen-
sion in order to find stochastic processes that are explicitly
solvable. We conclude in Section V with an outlook on fur-
ther possible developments.

## II. THE SHORT-RANGE EXPANSION OF THE FOKKER-PLANCK EQUATION WITH HETEROGENEOUS DIFFUSION

In this Section, will develop the formalism to study the short-time behavior of the one-dimensional Fokker-Planck equation. We initially state the problem under investigation and then we subdivide the development into two steps : (i) firstly, we obtain the singular short-time part of the propagator, and (ii) secondly we get the differential equations for the regular larger-time part. Finally, we obtain the formal solution of the hierarchy of differential equations describing the Taylor coefficients of the regular part of the propagator, particularly useful for the applications.

### A. Mathematical formulation of the problem

The starting point of our study is the Langevin equation, a stochastic differential equation of the form

$$ \frac{dx}{dt} = h(x,t) + g(x,t)\xi(t), \quad (1) $$

where $x(t)$ is the stochastic process under investigation, and $h$ and $g$ are regular functions. Without limiting generality, we can always assume that $g(x,t) \ge 0$ for any $x \in \mathbb{R}$ and $t \ge t_0$.

Moreover, $\xi(t)$ in Eq. (1) is a Gaussian white noise with average value $\mathbb{E}(\xi(t)) = 0$, and correlation $\mathbb{E}(\xi(t)\xi(t')) = 2\delta(t-t')$, where $\delta(t)$ is the Dirac delta function. Throughout this work, we always consider processes defined on the entire real axis, without any reflecting or absorbing boundary conditions.

For the process $x(t)$, we can introduce the probability density $W(x,t)$. It means that for any real numbers $x$ and $y$ (such that $y \le x$) we determine the probability that $y \le x(t) \le x$ as $\mathbb{P}(y \le x(t) \le x) = \int_y^x W(z,t)dz$, or equivalently, $W(x,t) = \frac{\partial \mathbb{P}(x(t) \le x)}{\partial x}$.

In order to completely define the meaning of Eq. (1), we must introduce a parameter $\alpha \in [0,1]$, identifying the integral stochastic interpretation adopted [1–3, 40, 41]. Indeed, the value of $\alpha$ defines the position of the point at which we calculate any integrated function in the small intervals of the adopted Riemann sum. This integration rule includes $\alpha = 0$ (Itô integral) [42], $\alpha = \frac{1}{2}$ (Fisk-Stratonovich integral) [43, 44], and $\alpha = 1$ (Hänggi-Klimontovich integral) [45, 46] (see also Ref. [47] for some comparison of these different choices).

The probability density $W(x,t)$ is the solution of the Fokker-Planck equation [1–3, 48–50]

$$ \frac{\partial W}{\partial t} = -\frac{\partial}{\partial x} \left[ \left( h + 2\alpha g \frac{\partial g}{\partial x} \right) W \right] + \frac{\partial^2}{\partial x^2} (g^2 W), \quad (2) $$

which can be also written in the equivalent form

$$ \frac{\partial W}{\partial t} = \frac{\partial}{\partial x} \left[ -hW + g^{2\alpha} \frac{\partial}{\partial x} \left( g^{2(1-\alpha)} W \right) \right]. \quad (3) $$

The right hand side of this Fokker-Planck equation can be expressed more explicitly as linear combination of $W$, $\frac{\partial W}{\partial x}$, and $\frac{\partial^2 W}{\partial x^2}$. We obtain

$$ \frac{\partial W}{\partial t} = -fW - \ell \frac{\partial W}{\partial x} + g^2 \frac{\partial^2 W}{\partial x^2} \quad (4) $$

where we introduced the regular functions

$$ f = \frac{\partial h}{\partial x} + 2(\alpha - 1) \left( \frac{\partial g}{\partial x} \right)^2 + 2(\alpha - 1) g \frac{\partial^2 g}{\partial x^2}, \quad (5) $$

$$ \ell = h + 2(\alpha - 2)g \frac{\partial g}{\partial x}. \quad (6) $$

These expressions show that the Fokker-Planck equation, and thus the density $W$, are not affected by the sign of the function $g$, which we assumed as always positive.

To solve the Fokker-Planck equation, we have to specify the initial condition $W(x,t_0) = W_0(x) \forall x \in \mathbb{R}$, at a fixed initial time $t_0$. In particular, a special solution is given by the kernel or propagator, defined by the initial condition $W(x,t_0) = \delta(x-y)$, where $y$ is an arbitrary initial value. It means that the initial value is deterministically known. This solution is referred to as $W(x,t) = K(x,t;y,t_0)$ and it can be obtained by the partial differential problem

$$ \frac{\partial K}{\partial t} = -fK - \ell \frac{\partial K}{\partial x} + g^2 \frac{\partial^2 K}{\partial x^2} \quad (7) $$

$$ K(x, t_0; y, t_0) = \delta(x - y). \quad (8) $$

With a little abuse of language, we can say that $K$ is the probability density of observing the value $x$ at time $t$, given the value $y$ at time $t_0$.

The propagator is important since any initial condition, such as $W(x,t_0) = W_0(x) \forall x \in \mathbb{R}$, can be handled by the relation

$$ W(x,t) = \int_{-\infty}^{+\infty} K(x,t;y,t_0)W_0(y)dy, \quad (9) $$

as a consequence of the linearity of the Fokker-Planck equation.

The objective of this work is to find a good approximation of the propagator $K$ for short times, and this will be organized as follows. The short-time approximation obviously contains the singularity induced by the initial condition represented by the delta function $\delta(x-y)$. It means that we look for $K_0$ such that $K(x,t;y,t_0) \sim K_0(x,t;y,t_0)$. Successively, for larger values of time, we can introduce a correction function $F$ such that $K(x,t;y,t_0) = K_0(x,t;y,t_0)F(x,t;y,t_0)$. Since the effects of the initial singularity are all contained in $K_0$, it follows that $F$ is regular. We will obtain a partial differen-
tial equation describing the behavior of $F$. Moreover, based on the regularity of $F$, we will introduce the Taylor development $F(x,t;y,t_0) = \sum_{n=0}^{+\infty} F_n(x;y)(t-t_0)^n$ and we obtain the coefficient $F_{n+1}$ recursively, i.e. as function of preceding coefficients $F_k, \forall k$ such that $0 \le k \le n$. Further simplifications are introduced in the case where the functions $h$ and $g$ are time-independent.

## B. Formal development of the theory

We describe here in detail the formal development of the theory, which includes the determination of the singular short-time part of the propagator and its regular short-time part. This regular perturbation is developed in series of powers over time, and the coefficients are described by differential equations that are solved explicitly.

The singular short-time part of the propagator. In order to obtain the explicit form of $K_0$ associated to the Fokker-Planck equation stated in Eq. (4), we firstly introduce the basic stochastic differential equation

$$ \frac{dx}{dt} = a\xi(t), \qquad (10) $$

where $a \in \mathbb{R}$. This represents the simplest Ornstein-Uhlenbeck process [51, 52], and the corresponding Fokker-Planck equation is given by the diffusion or heat equation

$$ \frac{\partial K}{\partial t} = a^2 \frac{\partial^2 K}{\partial x^2}. \qquad (11) $$

The propagator for this equation is well-known and can be written as follows [1, 2]

$$ K(x,t;y,t_0) = \frac{1}{\sqrt{4\pi a^2(t-t_0)}} \exp \left[ -\frac{(x-y)^2}{4a^2(t-t_0)} \right], \quad (12) $$

when it is determined from the initial condition $K(x,t_0;y,t_0) = \delta(x-y)$. This kernel represents the exact solution and perfectly describes the singularity for $t \to t_0$.

In the heat equation just solved, $a$ is a constant. In contrast, the diffusive term in the more general Fokker-Planck equation given in Eq. (4) arbitrarily depends on space and time. Hence, the idea is to consider the following short-time behavior of Eq. (4)

$$ K_0(x,t;y,t_0) = \frac{C}{\sqrt{4\pi(t-t_0)}} \exp \left[ -\frac{(\omega(x) - \omega(y))^2}{4(t-t_0)} \right]. \quad (13) $$

where $\omega(x)$ is a regular function able to take into account the heterogeneous diffusion.

The introduction of a function $\omega$, for now arbitrary, is useful to take into account the possible non-Gaussian behavior induced by the arbitrary form of $g(x,t)$. The validity of this proposed form of the short-time propagator is confirmed a posteriori, by substitution.

The constant $C$ has been introduced to normalize the propagator and is determined as follows, by considering the following representation of the delta function

$$ \lim_{\epsilon \to 0} \frac{1}{\sqrt{\pi \epsilon}} \exp\left(-\frac{z^2}{\epsilon}\right) = \delta(z). \qquad (14) $$

To this aim, we study the expression $\delta(\omega(x) - \omega(y))$ for an arbitrary $\omega(x)$. Of course, if $x \ne y$ we have $\delta(\omega(x) - \omega(y)) = 0$. Hence, we need to evaluate the integral $\int_{-\infty}^{+\infty} \delta(\omega(x) - \omega(y))dx$. To do this, we use the change of variable $z = \omega(x)$, and we introduce the inverse function $x = \Omega(z)$. We suppose that the function $\omega$ is always increasing, as we will be able to confirm later on. By definition, $\omega(\Omega(z)) = z$ and $\Omega(\omega(x)) = x$. Then, we get $dz = \omega'(x)dx$ and $dx = \Omega'(z)dz$, with $\omega'(x) \ge 0$ and $\Omega'(z) \ge 0$. We can write

$$ \begin{aligned} \int_{-\infty}^{+\infty} \delta(\omega(x) - \omega(y))dx &= \int_{-\infty}^{+\infty} \delta(z - \omega(y))\Omega'(z)dz \\ &= \Omega'(\omega(y)) = \frac{1}{\omega'(y)}. \end{aligned} \qquad (15) $$

Indeed, from $\Omega(\omega(x)) = x$, we obtain by derivation $\omega'(x)\Omega'(\omega(x)) = 1$ and hence $\Omega'(\omega(x)) = 1/\omega'(x)$. We finally proved that

$$ \delta(\omega(x) - \omega(y)) = \frac{1}{\omega'(y)} \delta(x - y). \qquad (16) $$

This is useful to obtain $C = \omega'(y)$, yielding the normalized expression for the propagator

$$ K_0(x,t;y,t_0) = \frac{\omega'(y)}{\sqrt{4\pi(t-t_0)}} \exp \left[ -\frac{(\omega(x) - \omega(y))^2}{4(t-t_0)} \right]. \quad (17) $$

Indeed, by using the classical limiting expression in Eq. (14), we deduce that

$$ K_0(x,t;y,t_0) \underset{t \to t_0}{\longrightarrow} \omega'(y)\delta(\omega(x) - \omega(y)) = \delta(x-y), \quad (18) $$

confirming the correct short time behavior of Eq. (17). The function $\omega(x)$ is useful to accommodate the complexity introduced by the space and time dependence of $f, \ell$, and in particular $g$. To obtain the explicit time-space dependence, we firstly substitute Eq. (17) into the Fokker-Planck equation. Hence, we need the following derivatives:

$$ \frac{\partial K_0}{\partial x} = -K_0 \frac{\omega(x) - \omega(y)}{2(t-t_0)} \omega'(x), \qquad (19) $$

$$ \frac{\partial^2 K_0}{\partial x^2} = K_0 \left[ \frac{\omega(x) - \omega(y)}{2(t-t_0)} \omega'(x) \right]^2 - \frac{K_0}{2(t-t_0)} \left[ \omega'(x)^2 + (\omega(x) - \omega(y))\omega''(x) \right], \quad (20) $$

$$ \frac{\partial K_0}{\partial t} = -\frac{K_0}{2(t-t_0)} + K_0 \frac{(\omega(x) - \omega(y))^2}{4(t-t_0)^2}. \qquad (21) $$

In Eq. (4), the leading terms for $t \to t_0$ are those proportional
tial equation describing the behavior of $F$. Moreover, based on the regularity of $F$, we will introduce the Taylor development $F(x,t;y,t_0) = \sum_{n=0}^{+\infty} F_n(x;y)(t-t_0)^n$ and we obtain the coefficient $F_{n+1}$ recursively, i.e. as function of preceding coefficients $F_k, \forall k$ such that $0 \le k \le n$. Further simplifications are introduced in the case where the functions $h$ and $g$ are time-independent.

## B. Formal development of the theory

We describe here in detail the formal development of the theory, which includes the determination of the singular short-time part of the propagator and its regular short-time part. This regular perturbation is developed in series of powers over time, and the coefficients are described by differential equations that are solved explicitly.

The singular short-time part of the propagator. In order to obtain the explicit form of $K_0$ associated to the Fokker-Planck equation stated in Eq. (4), we firstly introduce the basic stochastic differential equation

$$ \frac{dx}{dt} = a\xi(t), \qquad (10) $$

where $a \in \mathbb{R}$. This represents the simplest Ornstein-Uhlenbeck process [51, 52], and the corresponding Fokker-Planck equation is given by the diffusion or heat equation

$$ \frac{\partial K}{\partial t} = a^2 \frac{\partial^2 K}{\partial x^2}. \qquad (11) $$

The propagator for this equation is well-known and can be written as follows [1, 2]

$$ K(x,t;y,t_0) = \frac{1}{\sqrt{4\pi a^2(t-t_0)}} \exp \left[ -\frac{(x-y)^2}{4a^2(t-t_0)} \right], \quad (12) $$

when it is determined from the initial condition $K(x,t_0;y,t_0) = \delta(x-y)$. This kernel represents the exact solution and perfectly describes the singularity for $t \to t_0$.

In the heat equation just solved, $a$ is a constant. In contrast, the diffusive term in the more general Fokker-Planck equation given in Eq. (4) arbitrarily depends on space and time. Hence, the idea is to consider the following short-time behavior of Eq. (4)

$$ K_0(x,t;y,t_0) = \frac{C}{\sqrt{4\pi(t-t_0)}} \exp \left[ -\frac{(\omega(x) - \omega(y))^2}{4(t-t_0)} \right]. \quad (13) $$

where $\omega(x)$ is a regular function able to take into account the heterogeneous diffusion.

The introduction of a function $\omega$, for now arbitrary, is useful to take into account the possible non-Gaussian behavior induced by the arbitrary form of $g(x,t)$. The validity of this proposed form of the short-time propagator is confirmed a posteriori, by substitution.

The constant $C$ has been introduced to normalize the propagator and is determined as follows, by considering the following representation of the delta function

$$ \lim_{\epsilon \to 0} \frac{1}{\sqrt{\pi \epsilon}} \exp\left(-\frac{z^2}{\epsilon}\right) = \delta(z). \qquad (14) $$

To this aim, we study the expression $\delta(\omega(x) - \omega(y))$ for an arbitrary $\omega(x)$. Of course, if $x \ne y$ we have $\delta(\omega(x) - \omega(y)) = 0$. Hence, we need to evaluate the integral $\int_{-\infty}^{+\infty} \delta(\omega(x) - \omega(y))dx$. To do this, we use the change of variable $z = \omega(x)$, and we introduce the inverse function $x = \Omega(z)$. We suppose that the function $\omega$ is always increasing, as we will be able to confirm later on. By definition, $\omega(\Omega(z)) = z$ and $\Omega(\omega(x)) = x$. Then, we get $dz = \omega'(x)dx$ and $dx = \Omega'(z)dz$, with $\omega'(x) \ge 0$ and $\Omega'(z) \ge 0$. We can write

$$ \begin{aligned} \int_{-\infty}^{+\infty} \delta(\omega(x) - \omega(y))dx &= \int_{-\infty}^{+\infty} \delta(z - \omega(y))\Omega'(z)dz \\ &= \Omega'(\omega(y)) = \frac{1}{\omega'(y)}. \end{aligned} \qquad (15) $$

Indeed, from $\Omega(\omega(x)) = x$, we obtain by derivation $\omega'(x)\Omega'(\omega(x)) = 1$ and hence $\Omega'(\omega(x)) = 1/\omega'(x)$. We finally proved that

$$ \delta(\omega(x) - \omega(y)) = \frac{1}{\omega'(y)} \delta(x-y). \qquad (16) $$

This is useful to obtain $C = \omega'(y)$, yielding the normalized expression for the propagator

$$ K_0(x,t;y,t_0) = \frac{\omega'(y)}{\sqrt{4\pi(t-t_0)}} \exp \left[ -\frac{(\omega(x) - \omega(y))^2}{4(t-t_0)} \right]. \quad (17) $$

Indeed, by using the classical limiting expression in Eq. (14), we deduce that

$$ K_0(x,t;y,t_0) \underset{t \to t_0}{\longrightarrow} \omega'(y)\delta(\omega(x) - \omega(y)) = \delta(x-y), \quad (18) $$

confirming the correct short time behavior of Eq. (17). The function $\omega(x)$ is useful to accommodate the complexity introduced by the space and time dependence of $f, \ell$, and in particular $g$. To obtain the explicit time-space dependence, we firstly substitute Eq. (17) into the Fokker-Planck equation. Hence, we need the following derivatives:

$$ \frac{\partial K_0}{\partial x} = -K_0 \frac{\omega(x) - \omega(y)}{2(t-t_0)} \omega'(x), \qquad (19) $$

$$ \frac{\partial^2 K_0}{\partial x^2} = K_0 \left[ \frac{\omega(x) - \omega(y)}{2(t-t_0)} \omega'(x) \right]^2 - \frac{K_0}{2(t-t_0)} \left[ \omega'(x)^2 + (\omega(x) - \omega(y))\omega''(x) \right], \quad (20) $$

$$ \frac{\partial K_0}{\partial t} = -\frac{K_0}{2(t-t_0)} + K_0 \frac{(\omega(x) - \omega(y))^2}{4(t-t_0)^2}. \qquad (21) $$

In Eq. (4), the leading terms for $t \to t_0$ are those proportional
to $1/(t-t_0)^2$. We see that there is a leading term in $\frac{\partial^2 K_0}{\partial x^2}$, and another one in $\frac{\partial K_0}{\partial t}$. The balance of these two terms leads to

$$K_0 \frac{(\omega(x) - \omega(y))^2}{4(t-t_0)^2} = g(x, t_0)^2 K_0 \left[ \frac{\omega(x) - \omega(y)}{2(t-t_0)} \omega'(x) \right]^2, \quad (22)$$

and then $\omega'(x)^2 g^2(x, t_0) = 1$. Since $g(x, t_0) \ge 0$, we find $\omega'(x) = 1/g(x, t_0) \ge 0$, and therefore

$$\omega(x) - \omega(y) = \int_y^x \frac{d\eta}{g(\eta, t_0)} \triangleq \mathcal{D}(x, y, t_0), \quad (23)$$

where we defined the quantity $\mathcal{D}(x, y, t_0)$ for use in following developments. This calculation confirms that $\omega'(x) \ge 0$, as assumed in previous developments. To conclude, the short time propagator is given by

$$\begin{aligned} K_0(x,t;y,t_0) &= \frac{\exp\left[-\frac{1}{4(t-t_0)}\left(\int_y^x \frac{d\eta}{g(\eta,t_0)}\right)^2\right]}{g(y,t_0)\sqrt{4\pi(t-t_0)}} \\ &= \frac{\exp\left[-\frac{1}{4(t-t_0)}\mathcal{D}^2(x,y,t_0)\right]}{g(y,t_0)\sqrt{4\pi(t-t_0)}}, \end{aligned} \quad (24)$$

and it is only influenced by the diffusion coefficient $g$ for the initial time $t = t_0$. It is important to remark that the integral in Eq. (23) corresponds to the geodetic distance in a multi-dimensional setting, as discussed in Ref. [25]. Interestingly, the geodetic distance $\mathcal{D}$ can be written in closed form, as in Eq. (23), only for the one-dimensional case considered here.

The regular short-time part of the propagator. In order to improve the short-time representation of the propagator, we can consider the following solution

$$K(x,t;y,t_0) = K_0(x,t;y,t_0)F(x,t;y,t_0), \quad (25)$$

where $F$ is for the moment unknown. Of course, the choice to consider the multiplicative correction is arbitrary, and an additive correction could have been used. We kept the first choice because it involves analytically feasible calculations and exact solutions as discussed below and in the recent literature [23–25].

An equation for the quantity $F$ can be found by substituting Eq. (25) into Eq. (4), and by considering the expression just obtained for $K_0$. First of all, we introduce the notations $g_0(x) = g(x, t_0)$ and $g'_0(x) = \partial g(x, t_0)/\partial x$, and we calculate the following partial derivatives

$$\frac{\partial K}{\partial x} = -\frac{K_0}{2(t-t_0)} \frac{\mathcal{D}}{g_0} F + K_0 \frac{\partial F}{\partial x} \quad (26)$$

$$\frac{\partial K}{\partial t} = -\frac{K_0 F}{2(t-t_0)} + \frac{\mathcal{D}^2 K_0 F}{4(t-t_0)^2} + K_0 \frac{\partial F}{\partial t} \quad (27)$$

$$\begin{aligned} \frac{\partial^2 K}{\partial x^2} = & \frac{K_0}{4(t-t_0)^2} \frac{\mathcal{D}^2}{g_0^2} F - \frac{K_0}{2(t-t_0)} \frac{F}{g_0^2} + \frac{K_0}{2(t-t_0)} \frac{\mathcal{D} g'_0}{g_0^2} F \\ & - \frac{K_0}{(t-t_0)} \frac{\mathcal{D} \partial F}{g_0 \partial x} + K_0 \frac{\partial^2 F}{\partial x^2}. \end{aligned} \quad (28)$$

By using these expressions in Eq. (4), we elaborate an equation for the function $F$

$$\begin{aligned} \frac{\partial F}{\partial t} = & -fF - \ell \frac{\partial F}{\partial x} + g^2 \frac{\partial^2 F}{\partial x^2} - \frac{1}{2(t-t_0)} F \left( \frac{g_0^2}{g_0^2} - 1 \right) \\ & + \frac{1}{4(t-t_0)^2} \mathcal{D}^2 F \left( \frac{g_0^2}{g_0^2} - 1 \right) \\ & + \frac{\mathcal{D}}{2(t-t_0)g_0^2} \left( \ell F g_0 + g^2 g'_0 F - 2g^2 g_0 \frac{\partial F}{\partial x} \right). \end{aligned} \quad (29)$$

Although this result appears in a more complex form than the initial Fokker-Planck equation, the function $F$ always has a regular behavior for each time $t \ge t_0$. In such an equation, it must be taken care that $g(x, t)$ is different from $g_0(x) = g(x, t_0)$. We underline two specific points :

(i) there is only one initial condition to fully determine $F$, which is

$$F(y,t_0;y,t_0) = 1. \quad (30)$$

It means that the function $F$ must not alter the unit integral of the delta function for $t \to t_0$ and, therefore, it must be equal to one only for $x = y$.

(ii) since the singular behavior for $t \to t_0$ is fully described by $K_0$, the function $F$ is regular everywhere and can be developed in a Taylor series with respect to $t - t_0$.

According to this second point, since all the functions $F$, $f$, $\ell$, and $g^2$, involved in Eq. (29) are regular, we can introduce the following analytic developments

$$F(x,t;y,t_0) = \sum_{n=0}^{+\infty} F_n(x;y)(t-t_0)^n, \quad (31)$$

$$f(x,t) = \sum_{n=0}^{+\infty} f_n(x)(t-t_0)^n, \quad (32)$$

$$\ell(x,t) = \sum_{n=0}^{+\infty} \ell_n(x)(t-t_0)^n, \quad (33)$$

$$g^2(x,t) = \sum_{n=0}^{+\infty} g_n^2(x)(t-t_0)^n. \quad (34)$$

The initial condition stated in Eq. (30) simplifies to

$$F_0(y;y) = 1. \quad (35)$$

We will discuss later the conditions for $F_n$, $n \ge 1$, which need to be automatically identified by the procedure we are developing since there are no other a priori conditions to be imposed. We suppose from now on that all terms $f_n$, $\ell_n$, and $g_n^2$ are known regular functions of $x \in \mathbb{R}$, $n \ge 0$.

Before starting the substitution of Eqs. (31), (32), (33), and (34) into Eq. (29), we observe that

$$\frac{g^2}{g_0^2} - 1 = \sum_{n=1}^{+\infty} \frac{g_n^2}{g_0^2} (t - t_0)^n, \quad (36)$$

which means that $g^2/g_0^2 - 1$ is divisible by $t - t_0$. Therefore,
to $1/(t-t_0)^2$. We see that there is a leading term in $\frac{\partial^2 K_0}{\partial x^2}$, and another one in $\frac{\partial K_0}{\partial t}$. The balance of these two terms leads to

$$K_0 \frac{(\omega(x) - \omega(y))^2}{4(t - t_0)^2} = g(x, t_0)^2 K_0 \left[ \frac{\omega(x) - \omega(y)}{2(t - t_0)} \omega'(x) \right]^2, \quad (22)$$

and then $\omega'(x)^2 g^2(x, t_0) = 1$. Since $g(x, t_0) \ge 0$, we find $\omega'(x) = 1/g(x, t_0) \ge 0$, and therefore

$$\omega(x) - \omega(y) = \int_y^x \frac{d\eta}{g(\eta, t_0)} \triangleq \mathcal{D}(x, y, t_0), \quad (23)$$

where we defined the quantity $\mathcal{D}(x, y, t_0)$ for use in following developments. This calculation confirms that $\omega'(x) \ge 0$, as assumed in previous developments. To conclude, the short time propagator is given by

$$\begin{align}
K_0(x, t; y, t_0) &= \frac{\exp\left[-\frac{1}{4(t-t_0)}\left(\int_y^x \frac{d\eta}{g(\eta, t_0)}\right)^2\right]}{g(y, t_0)\sqrt{4\pi(t-t_0)}} \nonumber \\
&= \frac{\exp\left[-\frac{1}{4(t-t_0)}\mathcal{D}^2(x, y, t_0)\right]}{g(y, t_0)\sqrt{4\pi(t-t_0)}}, \tag{24}
\end{align}$$

and it is only influenced by the diffusion coefficient $g$ for the initial time $t = t_0$. It is important to remark that the integral in Eq. (23) corresponds to the geodetic distance in a multi-dimensional setting, as discussed in Ref. [25]. Interestingly, the geodetic distance $\mathcal{D}$ can be written in closed form, as in Eq. (23), only for the one-dimensional case considered here.

The regular short-time part of the propagator. In order to improve the short-time representation of the propagator, we can consider the following solution

$$K(x, t; y, t_0) = K_0(x, t; y, t_0) F(x, t; y, t_0), \quad (25)$$

where $F$ is for the moment unknown. Of course, the choice to consider the multiplicative correction is arbitrary, and an additive correction could have been used. We kept the first choice because it involves analytically feasible calculations and exact solutions as discussed below and in the recent literature [23–25].

An equation for the quantity $F$ can be found by substituting Eq. (25) into Eq. (4), and by considering the expression just obtained for $K_0$. First of all, we introduce the notations $g_0(x) = g(x, t_0)$ and $g'_0(x) = \partial g(x, t_0)/\partial x$, and we calculate the following partial derivatives

$$\frac{\partial K}{\partial x} = -\frac{K_0}{2(t-t_0)} \frac{\mathcal{D}}{g_0} F + K_0 \frac{\partial F}{\partial x} \quad (26)$$

$$\frac{\partial K}{\partial t} = -\frac{K_0 F}{2(t-t_0)} + \frac{\mathcal{D}^2 K_0 F}{4(t-t_0)^2} + K_0 \frac{\partial F}{\partial t} \quad (27)$$

$$\begin{align}
\frac{\partial^2 K}{\partial x^2} = & \frac{K_0}{4(t-t_0)^2} \frac{\mathcal{D}^2}{g_0^2} F - \frac{K_0}{2(t-t_0)} \frac{F}{g_0^2} + \frac{K_0}{2(t-t_0)} \frac{\mathcal{D} g'_0}{g_0^2} F \nonumber \\
& - \frac{K_0}{(t-t_0)} \frac{\mathcal{D} \partial F}{g_0 \partial x} + K_0 \frac{\partial^2 F}{\partial x^2}. \tag{28}
\end{align}$$

By using these expressions in Eq. (4), we elaborate an equation for the function $F$

$$\begin{align}
\frac{\partial F}{\partial t} = & -fF - \ell \frac{\partial F}{\partial x} + g^2 \frac{\partial^2 F}{\partial x^2} - \frac{1}{2(t-t_0)} F \left( \frac{g_0^2}{g^2} - 1 \right) \nonumber \\
& + \frac{1}{4(t-t_0)^2} \mathcal{D}^2 F \left( \frac{g^2}{g_0^2} - 1 \right) \tag{29} \\
& + \frac{\mathcal{D}}{2(t-t_0)g_0^2} \left( \ell F g_0 + g^2 g'_0 F - 2g^2 g_0 \frac{\partial F}{\partial x} \right). \nonumber
\end{align}$$

Although this result appears in a more complex form than the initial Fokker-Planck equation, the function $F$ always has a regular behavior for each time $t \ge t_0$. In such an equation, it must be taken care that $g(x, t)$ is different from $g_0(x) = g(x, t_0)$. We underline two specific points :

(i) there is only one initial condition to fully determine $F$, which is

$$F(y, t_0; y, t_0) = 1. \quad (30)$$

It means that the function $F$ must not alter the unit integral of the delta function for $t \to t_0$ and, therefore, it must be equal to one only for $x = y$.

(ii) since the singular behavior for $t \to t_0$ is fully described by $K_0$, the function $F$ is regular everywhere and can be developed in a Taylor series with respect to $t - t_0$.

According to this second point, since all the functions $F$, $f$, $\ell$, and $g^2$, involved in Eq. (29) are regular, we can introduce the following analytic developments

$$F(x, t; y, t_0) = \sum_{n=0}^{+\infty} F_n(x; y)(t - t_0)^n, \quad (31)$$

$$f(x, t) = \sum_{n=0}^{+\infty} f_n(x)(t - t_0)^n, \quad (32)$$

$$\ell(x, t) = \sum_{n=0}^{+\infty} \ell_n(x)(t - t_0)^n, \quad (33)$$

$$g^2(x, t) = \sum_{n=0}^{+\infty} g_n^2(x)(t - t_0)^n. \quad (34)$$

The initial condition stated in Eq. (30) simplifies to

$$F_0(y; y) = 1. \quad (35)$$

We will discuss later the conditions for $F_n$, $n \ge 1$, which need to be automatically identified by the procedure we are developing since there are no other a priori conditions to be imposed. We suppose from now on that all terms $f_n$, $\ell_n$, and $g_n^2$ are known regular functions of $x \in \mathbb{R}$, $n \ge 0$.

Before starting the substitution of Eqs. (31), (32), (33), and (34) into Eq. (29), we observe that

$$\frac{g^2}{g_0^2} - 1 = \sum_{n=1}^{+\infty} \frac{g_n^2}{g_0^2} (t - t_0)^n, \quad (36)$$

which means that $g^2/g_0^2 - 1$ is divisible by $t - t_0$. Therefore,
the last term in the first line of Eq. (29) is regular, i.e. without the singularity $1/(t-t_0)$, and the term in the second line exhibits a singularity of the type $1/(t-t_0)$ instead of type $1/(t-t_0)^2$, as it might seem at a first glance. Thus, the term in the second line together with the terms in the third line share the same singularity of the type $1/(t-t_0)$. Importantly, we will see that the balance equation for these terms will yield the first coefficient $F_0(x;y)$. That said, we can substitute all Taylor developments into the main equation and obtain, after long calculations, the relation

$$ \sum_{n=0}^{+\infty} F_{n+1}(n+1)(t-t_0)^n = \sum_{n=0}^{+\infty} \mathcal{G}_n(t-t_0)^n + \frac{\mathcal{D}^2 g_1^2}{4(t-t_0)g_0^2} F_0 + \frac{\mathcal{D} \ell_0}{2(t-t_0)g_0} F_0 + \frac{\mathcal{D} g_0'}{2(t-t_0)} F_0 - \frac{\mathcal{D} g_0}{(t-t_0)} \frac{\partial F_0}{\partial x}, \quad (37) $$

where

$$ \mathcal{G}_n = \sum_{k=0}^{n} \left[ -f_k F_{n-k} - \ell_k \frac{\partial F_{n-k}}{\partial x} + g_k^2 \frac{\partial^2 F_{n-k}}{\partial x^2} - \frac{g_{k+1}^2}{2g_0^2} F_{n-k} \right] + \sum_{k=0}^{n+1} \left[ \frac{\mathcal{D}^2}{4} \frac{g_{k+1}^2}{g_0^2} F_{n+1-k} + \frac{\mathcal{D}}{2g_0} \ell_k F_{n+1-k} + \frac{\mathcal{D} g_0'}{2g_0^2} g_k^2 F_{n+1-k} - \frac{\mathcal{D}}{g_0} g_k^2 \frac{\partial F_{n+1-k}}{\partial x} \right]. \quad (38) $$

The last four terms of Eq. (37), with the singularity $1/(t-t_0)$, give the following equation for $F_0$

$$ \frac{\partial F_0}{\partial x} = \left( \frac{\mathcal{D} g_1^2}{4 g_0^3} + \frac{\ell_0}{2 g_0^2} + \frac{g_0'}{2 g_0} \right) F_0, \quad (39) $$

where $F_0 = F_0(x;y)$, with $F_0(y;y) = 1$. Hence, the solution is easily obtained as

$$ F_0(x;y) = \exp \left[ \frac{1}{4} \int_y^x \frac{\mathcal{D}(\eta, y, t_0) g_1^2(\eta)}{g_0^3(\eta)} d\eta + \frac{1}{2} \int_y^x \frac{\ell_0(\eta)}{g_0^2(\eta)} d\eta + \frac{1}{2} \int_y^x \frac{g_0'(\eta)}{g_0(\eta)} d\eta \right]. \quad (40) $$

We therefore obtained a closed form expression for the first correction term $F_0$. The higher-order corrections can be found by equaling to zero the sum of terms having the same power $(t-t_0)^n$, $n \ge 0$, in Eq. (37). This results in a differential equation for $F_{n+1}$, $n \ge 0$, that depends on all the preceding terms

$F_0, \dots, F_n$. It reads

$$ \begin{aligned} & \mathcal{D} g_0 \frac{\partial F_{n+1}}{\partial x} + F_{n+1} \left( n+1 - \frac{\mathcal{D}^2 g_1^2}{4g_0^2} - \frac{\mathcal{D} \ell_0}{2g_0} - \frac{\mathcal{D} g_0'}{2} \right) \\ &= \sum_{k=0}^{n} \left( -f_k F_{n-k} - \ell_k \frac{\partial F_{n-k}}{\partial x} + g_k^2 \frac{\partial^2 F_{n-k}}{\partial x^2} - \frac{g_{k+1}^2}{2g_0^2} F_{n-k} \right. \\ & \qquad \left. + \frac{\mathcal{D}^2 g_{k+2}^2}{4g_0^2} F_{n-k} + \frac{\mathcal{D} \ell_{k+1}}{2g_0} F_{n-k} \right. \\ & \qquad \left. + \frac{\mathcal{D} g_0' g_{k+1}^2}{2g_0^2} F_{n-k} - \frac{\mathcal{D} g_{k+1}^2}{g_0} \frac{\partial F_{n-k}}{\partial x} \right). \end{aligned} \quad (41) $$

An important point concerns the initial condition to be used for this differential equation. The original condition $F(y,t_0;y,t_0) = 1$, simplified to $F_0(y;y) = 1$, has been already used in Eqs. (35) and (40). So, there are no other conditions from the original problem that we are studying. However, we easily see that the term with $\partial F_{n+1}/\partial x$ in Eq. (41) is proportional to $\mathcal{D}(x,y,t_0)$. It means that if we consider $x=y$ in Eq. (41) we obtain a relation fixing $F_{n+1}(y;y)$, since $\mathcal{D}(y,y,t_0) = 0$, see Eq. (23). The initial condition is therefore

$$ F_{n+1}(y;y) = \frac{1}{n+1} \sum_{k=0}^{n} \left( -f_k F_{n-k} - \ell_k \frac{\partial F_{n-k}}{\partial x} + g_k^2 \frac{\partial^2 F_{n-k}}{\partial x^2} - \frac{g_{k+1}^2}{2g_0^2} F_{n-k} \right). \quad (42) $$

From a conceptual point of view, the problem we studied was traced back to the solution of a hierarchy of differential problems, i.e. an infinite sequence of differential equations whose initial conditions are known. Nevertheless, although these equations are of the first order, their form is quite complex and needs to be simplified as described below.

### C. The formal solution of the hierarchy of differential equations

In order to simplify our results from section II B, we define

$$ \Phi(x,y) \equiv \frac{\mathcal{D} g_1^2}{4 g_0^3} + \frac{\ell_0}{2 g_0^2} + \frac{g_0'}{2 g_0}, \quad (43) $$

$$ \Psi_n(x,y) = \sum_{k=0}^{n-1} \left[ -f_k F_{n-1-k} - \ell_k \frac{\partial F_{n-1-k}}{\partial x} + g_k^2 \frac{\partial^2 F_{n-1-k}}{\partial x^2} - \frac{g_{k+1}^2}{2g_0^2} F_{n-1-k} + \frac{\mathcal{D}^2 g_{k+2}^2}{4g_0^2} F_{n-1-k} + \frac{\mathcal{D} \ell_{k+1}}{2g_0} F_{n-1-k} \right. \\ \left. + \frac{\mathcal{D} g_0' g_{k+1}^2}{2g_0^2} F_{n-1-k} - \frac{\mathcal{D} g_{k+1}^2}{g_0} \frac{\partial F_{n-1-k}}{\partial x} \right], \quad (44) $$
the last term in the first line of Eq. (29) is regular, i.e. without the singularity $1/(t-t_0)$, and the term in the second line exhibits a singularity of the type $1/(t-t_0)$ instead of type $1/(t-t_0)^2$, as it might seem at a first glance. Thus, the term in the second line together with the terms in the third line share the same singularity of the type $1/(t-t_0)$. Importantly, we will see that the balance equation for these terms will yield the first coefficient $F_0(x;y)$. That said, we can substitute all Taylor developments into the main equation and obtain, after long calculations, the relation

$$ \sum_{n=0}^{+\infty} F_{n+1}(n+1)(t-t_0)^n = \sum_{n=0}^{+\infty} \mathcal{G}_n(t-t_0)^n + \frac{\mathcal{D}^2 g_1^2}{4(t-t_0)g_0^2} F_0 + \frac{\mathcal{D} \ell_0}{2(t-t_0)g_0} F_0 + \frac{\mathcal{D} g_0'}{2(t-t_0)} F_0 - \frac{\mathcal{D} g_0}{(t-t_0)} \frac{\partial F_0}{\partial x}, \quad (37) $$

where

$$ \mathcal{G}_n = \sum_{k=0}^{n} \left[ -f_k F_{n-k} - \ell_k \frac{\partial F_{n-k}}{\partial x} + g_k^2 \frac{\partial^2 F_{n-k}}{\partial x^2} - \frac{g_{k+1}^2}{2g_0^2} F_{n-k} \right] + \sum_{k=0}^{n+1} \left[ \frac{\mathcal{D}^2 g_{k+1}^2}{4g_0^2} F_{n+1-k} + \frac{\mathcal{D} \ell_k}{2g_0} F_{n+1-k} + \frac{\mathcal{D} g_0'}{2g_0^2} g_k^2 F_{n+1-k} - \frac{\mathcal{D} g_0}{g_0} g_k^2 \frac{\partial F_{n+1-k}}{\partial x} \right]. \quad (38) $$

The last four terms of Eq. (37), with the singularity $1/(t-t_0)$, give the following equation for $F_0$

$$ \frac{\partial F_0}{\partial x} = \left( \frac{\mathcal{D} g_1^2}{4 g_0^3} + \frac{\ell_0}{2 g_0^2} + \frac{g_0'}{2 g_0} \right) F_0, \quad (39) $$

where $F_0 = F_0(x;y)$, with $F_0(y;y) = 1$. Hence, the solution is easily obtained as

$$ F_0(x;y) = \exp \left[ \frac{1}{4} \int_y^x \frac{\mathcal{D}(\eta, y, t_0) g_1^2(\eta)}{g_0^3(\eta)} d\eta + \frac{1}{2} \int_y^x \frac{\ell_0(\eta)}{g_0^2(\eta)} d\eta + \frac{1}{2} \int_y^x \frac{g_0'(\eta)}{g_0(\eta)} d\eta \right]. \quad (40) $$

We therefore obtained a closed form expression for the first correction term $F_0$. The higher-order corrections can be found by equaling to zero the sum of terms having the same power $(t-t_0)^n$, $n \ge 0$, in Eq. (37). This results in a differential equation for $F_{n+1}$, $n \ge 0$, that depends on all the preceding terms

$F_0, \dots, F_n$. It reads

$$ \begin{aligned} & \mathcal{D} g_0 \frac{\partial F_{n+1}}{\partial x} + F_{n+1} \left( n+1 - \frac{\mathcal{D}^2 g_1^2}{4g_0^2} - \frac{\mathcal{D} \ell_0}{2g_0} - \frac{\mathcal{D} g_0'}{2} \right) \\ &= \sum_{k=0}^{n} \left( -f_k F_{n-k} - \ell_k \frac{\partial F_{n-k}}{\partial x} + g_k^2 \frac{\partial^2 F_{n-k}}{\partial x^2} - \frac{g_{k+1}^2}{2g_0^2} F_{n-k} \right. \\ & \qquad \left. + \frac{\mathcal{D}^2 g_{k+2}^2}{4g_0^2} F_{n-k} + \frac{\mathcal{D} \ell_{k+1}}{2g_0} F_{n-k} \right. \\ & \qquad \left. + \frac{\mathcal{D} g_0' g_{k+1}^2}{2g_0^2} F_{n-k} - \frac{\mathcal{D} g_{k+1}^2}{g_0} \frac{\partial F_{n-k}}{\partial x} \right). \end{aligned} \quad (41) $$

An important point concerns the initial condition to be used for this differential equation. The original condition $F(y,t_0;y,t_0) = 1$, simplified to $F_0(y;y) = 1$, has been already used in Eqs. (35) and (40). So, there are no other conditions from the original problem that we are studying. However, we easily see that the term with $\partial F_{n+1}/\partial x$ in Eq. (41) is proportional to $\mathcal{D}(x,y,t_0)$. It means that if we consider $x=y$ in Eq. (41) we obtain a relation fixing $F_{n+1}(y;y)$, since $\mathcal{D}(y,y,t_0) = 0$, see Eq. (23). The initial condition is therefore

$$ F_{n+1}(y;y) = \frac{1}{n+1} \sum_{k=0}^{n} \left( -f_k F_{n-k} - \ell_k \frac{\partial F_{n-k}}{\partial x} + g_k^2 \frac{\partial^2 F_{n-k}}{\partial x^2} - \frac{g_{k+1}^2}{2g_0^2} F_{n-k} \right). \quad (42) $$

From a conceptual point of view, the problem we studied was traced back to the solution of a hierarchy of differential problems, i.e. an infinite sequence of differential equations whose initial conditions are known. Nevertheless, although these equations are of the first order, their form is quite complex and needs to be simplified as described below.

### C. The formal solution of the hierarchy of differential equations

In order to simplify our results from section II B, we define

$$ \Phi(x,y) \equiv \frac{\mathcal{D} g_1^2}{4 g_0^3} + \frac{\ell_0}{2 g_0^2} + \frac{g_0'}{2 g_0}, \quad (43) $$

$$ \Psi_n(x,y) = \sum_{k=0}^{n-1} \left[ -f_k F_{n-1-k} - \ell_k \frac{\partial F_{n-1-k}}{\partial x} + g_k^2 \frac{\partial^2 F_{n-1-k}}{\partial x^2} - \frac{g_{k+1}^2}{2g_0^2} F_{n-1-k} + \frac{\mathcal{D}^2 g_{k+2}^2}{4g_0^2} F_{n-1-k} + \frac{\mathcal{D} \ell_{k+1}}{2g_0} F_{n-1-k} \right. \\ \left. + \frac{\mathcal{D} g_0' g_{k+1}^2}{2g_0^2} F_{n-1-k} - \frac{\mathcal{D} g_{k+1}^2}{g_0} \frac{\partial F_{n-1-k}}{\partial x} \right], \quad (44) $$
and the differential equations, with corresponding initial conditions, can be rewritten as

$$ \frac{\partial F_0}{\partial x} - \Phi F_0 = 0, \quad (45) $$

$$ F_0(y, y) = 1, \quad (46) $$

$$ \frac{\partial F_n}{\partial x} + \left( n \frac{\mathcal{D}'}{\mathcal{D}} - \Phi \right) F_n = \frac{\mathcal{D}'}{\mathcal{D}} \Psi_n, \quad (47) $$

$$ F_n(y; y) = \frac{1}{n} \Psi_n(y, y), \quad (48) $$

where we considered $n \ge 1$, and we used the fact that $g_0(x) = 1/\mathcal{D}'(x,y,t_0)$, with $\mathcal{D}' = \partial\mathcal{D}/\partial x$. Since the quantity $\mathcal{D}$ depends on $y$, while $\mathcal{D}'$ does not depend on $y$, from now on, we use consistently the notation $\mathcal{D} = \mathcal{D}(x,y,t_0)$ and $\mathcal{D}' = \mathcal{D}'(x,t_0)$. Clearly, the solution of the differential equation for $F_0$ is

$$ F_0(x; y) = \exp \left[ \int_{y}^{x} \Phi(\eta, y) d\eta \right], \quad (49) $$

which exactly corresponds to Eq. (40). For solving Eq. (47), we start by considering the homogeneous equation $\partial F_n^h / \partial x + (n\mathcal{D}'/\mathcal{D} - \Phi) F_n^h = 0$, and we obtain the solution as

$$ \begin{aligned} F_n^h(x; y) &= \exp \left[ -n \int_z^x \frac{\mathcal{D}'(\eta, t_0)}{\mathcal{D}(\eta, y, t_0)} d\eta + \int_z^x \Phi(\eta, y) d\eta \right] \\ &= \exp \left[ -n \ln \left( \frac{\mathcal{D}(x, y, t_0)}{\mathcal{D}(z, y, t_0)} \right) + \int_z^x \Phi(\eta, y) d\eta \right] \\ &= \left( \frac{\mathcal{D}(z, y, t_0)}{\mathcal{D}(x, y, t_0)} \right)^n \exp \left[ \int_z^x \Phi(\eta, y) d\eta \right], \quad (50) \end{aligned} $$

where $\ln(\zeta)$ represents the natural logarithm of $\zeta$. Moreover, $z$ takes an arbitrary value that could be fixed by knowing the initial condition of $F_n^h$. For now, let's leave $z$ free and look for the solution of the non-homogeneous equation by adopting the method of variation of parameters (Lagrange's method). We can in fact write a particular solution in the form

$$ F_n^p(x; y) = C(x) \left( \frac{\mathcal{D}(z, y, t_0)}{\mathcal{D}(x, y, t_0)} \right)^n \exp \left[ \int_z^x \Phi(\eta, y) d\eta \right], \quad (51) $$

where $C(x)$ is an unknown regular function, which can be determined as follows. We can substitute Eq. (51) into Eq. (47) to get

$$ \begin{aligned} C'(x) \left( \frac{\mathcal{D}(z, y, t_0)}{\mathcal{D}(x, y, t_0)} \right)^n \exp \left[ \int_z^x \Phi(\eta, y) d\eta \right] \\ = \frac{\mathcal{D}'(x, t_0)}{\mathcal{D}(x, y, t_0)} \Psi_n(x, y). \end{aligned} \quad (52) $$

Hence, we find the following explicit expression for $C(x)$

$$ C(x) = \frac{1}{\mathcal{D}(z, y, t_0)^n} \int_y^x \mathcal{D}'(\eta, t_0) \mathcal{D}(\eta, y, t_0)^{n-1} \Psi_n(\eta, y) \times \exp \left[ -\int_z^\eta \Phi(\chi, y) d\chi \right] d\eta. \quad (53) $$

A particular solution of Eq. (47) is therefore

$$ \begin{aligned} F_n^p(x; y) = \frac{1}{\mathcal{D}(x, y, t_0)^n} \int_y^x \mathcal{D}'(\eta, t_0) \mathcal{D}(\eta, y, t_0)^{n-1} \Psi_n(\eta, y) \\ \times \exp \left[ -\int_z^\eta \Phi(\chi, y) d\chi \right] \exp \left[ \int_z^x \Phi(\eta, y) d\eta \right], \end{aligned} \quad (54) $$

which is valid for any value of $z$. In general, following Lagrange theory, the complete solution of Eq. (47) is given by the sum of the solution of the homogeneous equation and the particular solution. In this specific case, the procedure can be simplified by trying to consider $z = y$ in previous solutions. So doing, the solution of the homogeneous equation is zero since $\mathcal{D}(y, y, t_0) = 0$, and the general solution becomes coincident with the particular solution with $z = y$ ($n \ge 1$)

$$ F_n(x; y) = \frac{1}{\mathcal{D}(x, y, t_0)^n} \int_y^x \mathcal{D}'(\eta, t_0) \mathcal{D}(\eta, y, t_0)^{n-1} \times \Psi_n(\eta, y) \exp \left[ -\int_y^\eta \Phi(\chi, y) d\chi \right] \exp \left[ \int_z^x \Phi(\eta, y) d\eta \right]. \quad (55) $$

We can verify that this solution perfectly satisfies the initial condition $F_n(y; y) = \Psi_n(y, y)/n$. To this aim, we consider two regular functions $\varphi(x)$ and $a(x)$ such that $\varphi(y) = 0$, and we use the result

$$ \lim_{x \to y} \frac{1}{\varphi(x)} \int_y^x a(\eta) \varphi'(\eta) d\eta = a(y), \quad (56) $$

which can be directly proved, for instance, by means of L'Hôpital's rule. We can assume $\varphi(x) = \mathcal{D}(x, y, t_0)^n$ so that $\varphi'(x) = n\mathcal{D}'(x, t_0)\mathcal{D}(x, y, t_0)^{n-1}$, and $a(x) = \frac{1}{n}\Psi_n(x, y)\exp(-\int_y^x \Phi(\chi, y) d\chi)$. Then, we find that $\lim_{x \to y} F_n(x; y) = \Psi_n(y, y)/n$, as expected.

To conclude, Eqs. (49) and (55) are the formal solutions of the hierarchy of differential equations for the coefficients $F_n$. Thus we have at our disposal the complete mathematical form of the propagator $K = K_0 F$ of the Fokker-Planck equation, consisting in the product of Eq. (24) and Eq. (31), where the coefficients just found are to be used.

A particular but very important case concerns the time-independent stochastic differential equation

$$ \frac{dx}{dt} = h(x) + g(x) \xi(t), \quad (57) $$

In this case the functions $f$ and $\ell$ remains unchanged with respect to Eqs. (5) and (6), but their Taylor developments are greatly simplified. Indeed, we have $f_0(x) = f(x)$, $\ell_0(x) = \ell(x)$, and $g_0^2(x) = g^2(x)$. Moreover, for all $n \ge 1$, we have $f_n(x) = 0$, $\ell_n(x) = 0$, and $g_n^2(x) = 0$. As a result, the functions $\Phi$ and $\Psi_n$
and the differential equations, with corresponding initial conditions, can be rewritten as

$$ \frac{\partial F_0}{\partial x} - \Phi F_0 = 0, \quad (45) $$

$$ F_0(y, y) = 1, \quad (46) $$

$$ \frac{\partial F_n}{\partial x} + \left( n \frac{\mathcal{D}'}{\mathcal{D}} - \Phi \right) F_n = \frac{\mathcal{D}'}{\mathcal{D}} \Psi_n, \quad (47) $$

$$ F_n(y; y) = \frac{1}{n} \Psi_n(y, y), \quad (48) $$

where we considered $n \ge 1$, and we used the fact that $g_0(x) = 1/\mathcal{D}'(x,y,t_0)$, with $\mathcal{D}' = \partial\mathcal{D}/\partial x$. Since the quantity $\mathcal{D}$ depends on $y$, while $\mathcal{D}'$ does not depend on $y$, from now on, we use consistently the notation $\mathcal{D} = \mathcal{D}(x,y,t_0)$ and $\mathcal{D}' = \mathcal{D}'(x,t_0)$. Clearly, the solution of the differential equation for $F_0$ is

$$ F_0(x; y) = \exp \left[ \int_{y}^{x} \Phi(\eta, y) d\eta \right], \quad (49) $$

which exactly corresponds to Eq. (40). For solving Eq. (47), we start by considering the homogeneous equation $\partial F_n^h / \partial x + (n\mathcal{D}'/\mathcal{D} - \Phi) F_n^h = 0$, and we obtain the solution as

$$ \begin{aligned} F_n^h(x; y) &= \exp \left[ -n \int_z^x \frac{\mathcal{D}'(\eta, t_0)}{\mathcal{D}(\eta, y, t_0)} d\eta + \int_z^x \Phi(\eta, y) d\eta \right] \\ &= \exp \left[ -n \ln \left( \frac{\mathcal{D}(x, y, t_0)}{\mathcal{D}(z, y, t_0)} \right) + \int_z^x \Phi(\eta, y) d\eta \right] \\ &= \left( \frac{\mathcal{D}(z, y, t_0)}{\mathcal{D}(x, y, t_0)} \right)^n \exp \left[ \int_z^x \Phi(\eta, y) d\eta \right], \quad (50) \end{aligned} $$

where $\ln(\zeta)$ represents the natural logarithm of $\zeta$. Moreover, $z$ takes an arbitrary value that could be fixed by knowing the initial condition of $F_n^h$. For now, let's leave $z$ free and look for the solution of the non-homogeneous equation by adopting the method of variation of parameters (Lagrange's method). We can in fact write a particular solution in the form

$$ F_n^p(x; y) = C(x) \left( \frac{\mathcal{D}(z, y, t_0)}{\mathcal{D}(x, y, t_0)} \right)^n \exp \left[ \int_z^x \Phi(\eta, y) d\eta \right], \quad (51) $$

where $C(x)$ is an unknown regular function, which can be determined as follows. We can substitute Eq. (51) into Eq. (47) to get

$$ \begin{aligned} C'(x) \left( \frac{\mathcal{D}(z, y, t_0)}{\mathcal{D}(x, y, t_0)} \right)^n \exp \left[ \int_z^x \Phi(\eta, y) d\eta \right] \\ = \frac{\mathcal{D}'(x, t_0)}{\mathcal{D}(x, y, t_0)} \Psi_n(x, y). \end{aligned} \quad (52) $$

Hence, we find the following explicit expression for $C(x)$

$$ C(x) = \frac{1}{\mathcal{D}(z, y, t_0)^n} \int_y^x \mathcal{D}'(\eta, t_0) \mathcal{D}(\eta, y, t_0)^{n-1} \Psi_n(\eta, y) \times \exp \left[ -\int_z^\eta \Phi(\chi, y) d\chi \right] d\eta. \quad (53) $$

A particular solution of Eq. (47) is therefore

$$ \begin{aligned} F_n^p(x; y) = \frac{1}{\mathcal{D}(x, y, t_0)^n} \int_y^x \mathcal{D}'(\eta, t_0) \mathcal{D}(\eta, y, t_0)^{n-1} \Psi_n(\eta, y) \\ \times \exp \left[ -\int_z^\eta \Phi(\chi, y) d\chi \right] \exp \left[ \int_z^x \Phi(\eta, y) d\eta \right], \end{aligned} \quad (54) $$

which is valid for any value of $z$. In general, following Lagrange theory, the complete solution of Eq. (47) is given by the sum of the solution of the homogeneous equation and the particular solution. In this specific case, the procedure can be simplified by trying to consider $z = y$ in previous solutions. So doing, the solution of the homogeneous equation is zero since $\mathcal{D}(y, y, t_0) = 0$, and the general solution becomes coincident with the particular solution with $z = y$ ($n \ge 1$)

$$ F_n(x; y) = \frac{1}{\mathcal{D}(x, y, t_0)^n} \int_y^x \mathcal{D}'(\eta, t_0) \mathcal{D}(\eta, y, t_0)^{n-1} \times \Psi_n(\eta, y) \exp \left[ -\int_y^\eta \Phi(\chi, y) d\chi \right] \exp \left[ \int_z^x \Phi(\eta, y) d\eta \right]. \quad (55) $$

We can verify that this solution perfectly satisfies the initial condition $F_n(y; y) = \Psi_n(y, y)/n$. To this aim, we consider two regular functions $\varphi(x)$ and $a(x)$ such that $\varphi(y) = 0$, and we use the result

$$ \lim_{x \to y} \frac{1}{\varphi(x)} \int_y^x a(\eta) \varphi'(\eta) d\eta = a(y), \quad (56) $$

which can be directly proved, for instance, by means of L'Hôpital's rule. We can assume $\varphi(x) = \mathcal{D}(x, y, t_0)^n$ so that $\varphi'(x) = n\mathcal{D}'(x, t_0)\mathcal{D}(x, y, t_0)^{n-1}$, and $a(x) = \frac{1}{n}\Psi_n(x, y)\exp(-\int_y^x \Phi(\chi, y) d\chi)$. Then, we find that $\lim_{x \to y} F_n(x; y) = \Psi_n(y, y)/n$, as expected.

To conclude, Eqs. (49) and (55) are the formal solutions of the hierarchy of differential equations for the coefficients $F_n$. Thus we have at our disposal the complete mathematical form of the propagator $K = K_0 F$ of the Fokker-Planck equation, consisting in the product of Eq. (24) and Eq. (31), where the coefficients just found are to be used.

A particular but very important case concerns the time-independent stochastic differential equation

$$ \frac{dx}{dt} = h(x) + g(x) \xi(t), \quad (57) $$

In this case the functions $f$ and $\ell$ remains unchanged with respect to Eqs. (5) and (6), but their Taylor developments are greatly simplified. Indeed, we have $f_0(x) = f(x)$, $\ell_0(x) = \ell(x)$, and $g_0^2(x) = g^2(x)$. Moreover, for all $n \ge 1$, we have $f_n(x) = 0$, $\ell_n(x) = 0$, and $g_n^2(x) = 0$. As a result, the functions $\Phi$ and $\Psi_n$
turn out to be simplified as follows

$$ \Phi = \frac{\ell}{2g^2} + \frac{g'}{2g}, \qquad (58) $$

$$ \Psi_n = -f F_{n-1} - \ell \frac{\partial F_{n-1}}{\partial x} + g^2 \frac{\partial^2 F_{n-1}}{\partial x^2} \quad (59) $$

Importantly, $\Psi_n$ depends only on the previous coefficient $F_{n-1}$ and not on all coefficients as in the arbitrarily time-varying case. So, in this case Eq. (55) represents a simple recursion between $F_n$ and $F_{n-1}$. This makes the procedure easily applicable in many practical cases. In addition, this approach can be further simplified by defining

$$ F_n(x; y) = D_n(x, y) F_0(x; y). \qquad (60) $$

Of course, we have $D_0 = 1$. As discussed below, the recursion law for the parameters $D_n$ is simpler than the previous one, concerning the parameters $F_n$. To prove this point, we begin by developing the derivatives in Eq. (59), as follows

$$ \frac{\partial F_k}{\partial x} = \left( \frac{\partial D_k}{\partial x} + D_k \Phi \right) F_0 \qquad (61) $$

$$ \frac{\partial^2 F_k}{\partial x^2} = \left( \frac{\partial^2 D_k}{\partial x^2} + 2 \frac{\partial D_k}{\partial x} \Phi + D_k \frac{\partial \Phi}{\partial x} + D_k \Phi^2 \right) F_0. \quad (62) $$

It is important to remark that both results are proportional to $F_0$. It means that we can introduce a function $\Lambda_n$, such that $\Psi_n = \Lambda_n F_0$, which can be eventually written as

$$ \Lambda_n = -f D_{n-1} - \ell \left( \frac{\partial D_{n-1}}{\partial x} + D_{n-1} \Phi \right) + g^2 \left( \frac{\partial^2 D_{n-1}}{\partial x^2} + 2 \frac{\partial D_{n-1}}{\partial x} \Phi + D_{n-1} \frac{\partial \Phi}{\partial x} + D_{n-1} \Phi^2 \right). \quad (63) $$

In this time-independent situation, the function $\mathcal{D}$ can be redefined as $\mathcal{D}(x,y) = \int_y^x d\eta/g(\eta)$ and, therefore, $\mathcal{D}'(x) = 1/g(x)$. Hence, Eq. (55) takes the form

$$ \begin{aligned} F_n(x; y) &= \int_y^x \frac{\mathcal{D}'(\eta) \mathcal{D}(\eta, y)^{n-1} \Psi_n(\eta, y)}{\mathcal{D}(x, y)^n F_0(\eta, y)} d\eta F_0(x; y) \\ &= \int_y^x \frac{\mathcal{D}'(\eta) \mathcal{D}(\eta, y)^{n-1} \Lambda_n(\eta, y)}{\mathcal{D}(x, y)^n} d\eta F_0(x; y), \end{aligned} \quad (64) $$

or, equivalently

$$ D_n(x, y) = \frac{1}{\mathcal{D}(x, y)^n} \int_y^x \mathcal{D}'(\eta) \mathcal{D}(\eta, y)^{n-1} \Lambda_n(\eta, y) d\eta. \quad (65) $$

Finally, the combination of Eqs. (63) and (65) represents a recursive procedure that gives all the parameters $D_n$, which, in turn, can be used to obtain all the parameters $F_n$ through Eq. (60). This approach can be easily applied to study the behavior of stochastic differential equations, as in Eq. (57), with arbitrary drift and heterogeneous diffusion.

In the next section, we will turn to applications of our

formalism. In particular, we demonstrate the determination of the coefficients $D_n$, and the convergence of the series in Eq. (31) for various stochastic processes.

### III. APPLICATION TO SPECIFIC STOCHASTIC PROCESSES

We now turn to the application of our formal developments to specific examples. We begin with the simplest of all cases, a Gaussian process, which has numerous applications in many scientific fields. The second is the geometric Brownian process, a stochastic process most often associated with problems in finance. It is worth noting that for these first two examples we could compute the full propagator exactly, against which the form of short-term propagator can be easily compared. The third and fourth are two case studies taken from biophysics, for which we do not know a priori the full propagator: (i) a model stochastic process for a molecular motor, and (ii) a stochastic process occurring in parasite motility, based on an exponential heterogeneous diffusion coefficient.

#### A. Gaussian processes

This simple case describes a system with constant drift and diffusion terms. We consider therefore two real parameters $H_0$ and $G_0 > 0$ and the stochastic differential equation

$$ \frac{dx}{dt} = -H_0 + G_0 \xi(t). \qquad (66) $$

Of course, the propagator for this problem is well-known, in fact a generalization of Eq. 12 that can be written as [1, 2]

$$ K(x, t; y, t_0) = \frac{\exp \left[ - \frac{(x-y+H_0(t-t_0))^2}{4G_0^2(t-t_0)} \right]}{\sqrt{4\pi G_0^2(t-t_0)}}. \quad (67) $$

Within our formalism we thus have for Eqs. (5) and (6) $f(x) = 0$, $\ell(x) = -H_0$, and $g^2(x) = G_0^2$. Moreover, $\mathcal{D}(x, y) = (x-y)/G_0$. Therefore, we can obtain the short-time propagator $K_0$ in the form

$$ K_0(x, t; y, t_0) = \frac{1}{\sqrt{4\pi G_0^2(t - t_0)}} \exp \left[ - \frac{(x - y)^2}{4G_0^2(t - t_0)} \right]. \quad (68) $$

It can be immediately seen that this expression, valid for short times, does not coincide with $K$, Eq. (67). Therefore, the correction introduced by the function $F$ is essential in this case to obtain the correct propagator.

To determine this correction, we evaluate $\Phi = -\frac{H_0}{2G_0^2}$, and we write $F_0$ as

$$ F_0(x; y) = \exp \left[ \int_y^x \Phi(\eta) d\eta \right] = \exp \left[ - \frac{H_0}{2G_0^2} (x - y) \right]. \quad (69) $$
turn out to be simplified as follows

$$ \Phi = \frac{\ell}{2g^2} + \frac{g'}{2g}, \qquad (58) $$

$$ \Psi_n = -f F_{n-1} - \ell \frac{\partial F_{n-1}}{\partial x} + g^2 \frac{\partial^2 F_{n-1}}{\partial x^2} \quad (59) $$

Importantly, $\Psi_n$ depends only on the previous coefficient $F_{n-1}$ and not on all coefficients as in the arbitrarily time-varying case. So, in this case Eq. (55) represents a simple recursion between $F_n$ and $F_{n-1}$. This makes the procedure easily applicable in many practical cases. In addition, this approach can be further simplified by defining

$$ F_n(x; y) = D_n(x, y) F_0(x; y). \qquad (60) $$

Of course, we have $D_0 = 1$. As discussed below, the recursion law for the parameters $D_n$ is simpler than the previous one, concerning the parameters $F_n$. To prove this point, we begin by developing the derivatives in Eq. (59), as follows

$$ \frac{\partial F_k}{\partial x} = \left( \frac{\partial D_k}{\partial x} + D_k \Phi \right) F_0 \qquad (61) $$

$$ \frac{\partial^2 F_k}{\partial x^2} = \left( \frac{\partial^2 D_k}{\partial x^2} + 2 \frac{\partial D_k}{\partial x} \Phi + D_k \frac{\partial \Phi}{\partial x} + D_k \Phi^2 \right) F_0. \quad (62) $$

It is important to remark that both results are proportional to $F_0$. It means that we can introduce a function $\Lambda_n$, such that $\Psi_n = \Lambda_n F_0$, which can be eventually written as

$$ \Lambda_n = -f D_{n-1} - \ell \left( \frac{\partial D_{n-1}}{\partial x} + D_{n-1} \Phi \right) + g^2 \left( \frac{\partial^2 D_{n-1}}{\partial x^2} + 2 \frac{\partial D_{n-1}}{\partial x} \Phi + D_{n-1} \frac{\partial \Phi}{\partial x} + D_{n-1} \Phi^2 \right). \quad (63) $$

In this time-independent situation, the function $\mathcal{D}$ can be redefined as $\mathcal{D}(x,y) = \int_y^x d\eta/g(\eta)$ and, therefore, $\mathcal{D}'(x) = 1/g(x)$. Hence, Eq. (55) takes the form

$$ \begin{aligned} F_n(x; y) &= \int_y^x \frac{\mathcal{D}'(\eta) \mathcal{D}(\eta, y)^{n-1} \Psi_n(\eta, y)}{\mathcal{D}(x, y)^n F_0(\eta, y)} d\eta F_0(x; y) \\ &= \int_y^x \frac{\mathcal{D}'(\eta) \mathcal{D}(\eta, y)^{n-1} \Lambda_n(\eta, y)}{\mathcal{D}(x, y)^n} d\eta F_0(x; y), \end{aligned} \quad (64) $$

or, equivalently

$$ D_n(x, y) = \frac{1}{\mathcal{D}(x, y)^n} \int_y^x \mathcal{D}'(\eta) \mathcal{D}(\eta, y)^{n-1} \Lambda_n(\eta, y) d\eta. \quad (65) $$

Finally, the combination of Eqs. (63) and (65) represents a recursive procedure that gives all the parameters $D_n$, which, in turn, can be used to obtain all the parameters $F_n$ through Eq. (60). This approach can be easily applied to study the behavior of stochastic differential equations, as in Eq. (57), with arbitrary drift and heterogeneous diffusion.

In the next section, we will turn to applications of our

formalism. In particular, we demonstrate the determination of the coefficients $D_n$, and the convergence of the series in Eq. (31) for various stochastic processes.

### III. APPLICATION TO SPECIFIC STOCHASTIC PROCESSES

We now turn to the application of our formal developments to specific examples. We begin with the simplest of all cases, a Gaussian process, which has numerous applications in many scientific fields. The second is the geometric Brownian process, a stochastic process most often associated with problems in finance. It is worth noting that for these first two examples we could compute the full propagator exactly, against which the form of short-term propagator can be easily compared. The third and fourth are two case studies taken from biophysics, for which we do not know a priori the full propagator: (i) a model stochastic process for a molecular motor, and (ii) a stochastic process occurring in parasite motility, based on an exponential heterogeneous diffusion coefficient.

#### A. Gaussian processes

This simple case describes a system with constant drift and diffusion terms. We consider therefore two real parameters $H_0$ and $G_0 > 0$ and the stochastic differential equation

$$ \frac{dx}{dt} = -H_0 + G_0 \xi(t). \qquad (66) $$

Of course, the propagator for this problem is well-known, in fact a generalization of Eq. 12 that can be written as [1, 2]

$$ K(x, t; y, t_0) = \frac{\exp \left[ - \frac{(x-y+H_0(t-t_0))^2}{4G_0^2(t-t_0)} \right]}{\sqrt{4\pi G_0^2(t-t_0)}}. \quad (67) $$

Within our formalism we thus have for Eqs. (5) and (6) $f(x) = 0$, $\ell(x) = -H_0$, and $g^2(x) = G_0^2$. Moreover, $\mathcal{D}(x, y) = (x - y)/G_0$. Therefore, we can obtain the short-time propagator $K_0$ in the form

$$ K_0(x, t; y, t_0) = \frac{1}{\sqrt{4\pi G_0^2(t - t_0)}} \exp \left[ - \frac{(x - y)^2}{4G_0^2(t - t_0)} \right]. \quad (68) $$

It can be immediately seen that this expression, valid for short times, does not coincide with $K$, Eq. (67). Therefore, the correction introduced by the function $F$ is essential in this case to obtain the correct propagator.

To determine this correction, we evaluate $\Phi = -\frac{H_0}{2G_0^2}$, and we write $F_0$ as

$$ F_0(x; y) = \exp \left[ \int_y^x \Phi(\eta) d\eta \right] = \exp \left[ - \frac{H_0}{2G_0^2} (x - y) \right]. \quad (69) $$
<figure><img src="image_1.png" /><figcaption>FIGURE 1. Demonstration of the convergence of the expansion for a Gaussian process. We plot the short-time propagator $K_0$, the full propagator $K$, and low-order approximations $K_0Q_0$, $K_0Q_1$, $K_0Q_2$, and $K_0Q_3$, to show the convergence $K_0Q_m \to K$ when $m \to \infty$. We adopted the parameter values $H_0 = -1$, $G_0 = 1$, $y = 1$, $t_0 = 0$, $t = 3$ (in arbitrary units).</figcaption></figure>

We can prove by mathematical induction that

$$D_n = \frac{(-1)^n H_0^{2n}}{4^n n! G_0^{2n}}. \qquad (70)$$

This expression was found by explicitly applying the recursive procedure - see Eqs.(63) and (65) - to the first terms and looking for a general formula, which must now be rigorously proved. To begin with, we note that the formula gives the correct value $D_0 = 1$, for $n = 0$. Now, we can calculate $D_{n+1}$ with the procedure of Section IIC to verify that the assumption made on $D_n$ is correct. To this aim, we firstly determine $\Lambda_{n+1}$ through Eq. (63) as follows

$$\Lambda_{n+1} = \frac{(-1)^{n+1} H_0^{2(n+1)}}{4^{n+1} n! G_0^{2(n+1)}}. \qquad (71)$$

Then, we can determine $D_{n+1}$ by means of Eq. (65)

$$ \begin{aligned} D_{n+1} &= \int_y^x \frac{\mathcal{D}'(\eta)\mathcal{D}(\eta, y)^n}{\mathcal{D}(x, y)^{n+1}} \frac{(-1)^{n+1} H_0^{2(n+1)}}{4^{n+1} n! G_0^{2(n+1)}} d\eta \\ &= \int_y^x \frac{(\eta - y)^n}{(x - y)^{n+1}} \frac{(-1)^{n+1} H_0^{2(n+1)}}{4^{n+1} n! G_0^{2(n+1)}} d\eta \\ &= \frac{(-1)^{n+1} H_0^{2(n+1)}}{4^{n+1} (n+1)! G_0^{2(n+1)}}, \end{aligned} \qquad (72) $$

which is consistent with Eq. (70) through the replacement $n \to n+1$. This confirms the validity of the result by mathematical induction.

This approach, based on an *ansatz* like that in Eq. (70), which is then verified by induction on the natural numbers is interesting but can probably be applied to a variety of special cases but not to the general problem. In fact, in the most com-

plicated cases it is difficult to propose an expression for $D_n$ by applying only a few recursions of the procedure. However, it remains useful for various applications discussed below.

We can now explicitly determine $F$ as follows

$$ \begin{aligned} F(x, t; y, t_0) &= \sum_{n=0}^{+\infty} F_n(x; y)(t - t_0)^n \\ &= \sum_{n=0}^{+\infty} D_n(x)(t - t_0)^n F_0(x; y) \\ &= \sum_{n=0}^{+\infty} \frac{\left(-\frac{H_0^2(t-t_0)}{4G_0^2}\right)^n}{n!} F_0(x; y) \\ &= \exp\left[-\frac{H_0^2(t-t_0)}{4G_0^2}\right] \exp\left[-\frac{H_0}{2G_0^2}(x-y)\right]. \end{aligned} \qquad (73) $$

We can finally multiply $K_0$ by $F$, and it is easily seen that we get the correct solution for $K$, given by Eq. (67). The series defining the function $F$ is clearly convergent and thus determines the first exponential function in Eq. (73). Although rather simple, this example shows that the coefficients $D_n$ decrease rapidly with $n$ and thus a few terms are often sufficient to have a good approximation on short-time intervals.

An illustration of this behaviour is found in Fig. 1, where we show $K_0$, $K$, and different approximations given by $K_0Q_m$, where $Q_m = \sum_{n=0}^m F_n(x; y)(t - t_0)^n$; $Q_m$ is partial sum of Eq. (73). We see that the first approximation $K_0Q_0$ is too crude and in fact is not even properly normalizable. On the other hand, the subsequent approximations are automatically normalized to within a negligible error. To conclude, we finally remark that the results of this first example are independent of $\alpha$ since the noise is additive and not multiplicative. We next turn to a first example with multiplicative noise.

## B. Geometric Brownian processes

Geometric Brownian processes find application in physics, e.g., in the statistical theory of turbulence [38, 39], but they are much more popular in finance, in particular in the modeling of stock option prices [53–56].

For this example we will see that we can compute the kernel with our recursive procedure for an arbitrary interpretation of the stochastic integrals (i.e. for any value of $\alpha$). We consider an exemplary geometric Brownian process described by the stochastic differential equation

$$\frac{dx}{dt} = -H_0x + G_0|x|\xi(t), \qquad (74)$$

where we considered the absolute value in the diffusion term in order to have $g(x) \ge 0$, $x \in \mathbb{R}$, that is a key assumption of all our work (considering $G_0 > 0$). We have $h(x) = -H_0x$, $g(x) = G_0|x|$, and $g'(x) = G_0 \operatorname{sgn}(x)$, where $\operatorname{sgn}(\zeta)$ represents the signum function returning the sign of $\zeta$, i.e., $\operatorname{sgn}(\zeta) = +1$ if $\zeta > 0$, $\operatorname{sgn}(\zeta) = 0$ if $\zeta = 0$, and $\operatorname{sgn}(\zeta) = -1$ if $\zeta < 0$. It
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Demonstration of the convergence of the expansion for a Gaussian process. We plot the short-time propagator $K_0$, the full propagator $K$, and low-order approximations $K_0Q_0$, $K_0Q_1$, $K_0Q_2$, and $K_0Q_3$, to show the convergence $K_0Q_m \to K$ when $m \to \infty$. We adopted the parameter values $H_0 = -1$, $G_0 = 1$, $y = 1$, $t_0 = 0$, $t = 3$ (in arbitrary units).</figcaption></figure>

We can prove by mathematical induction that

$$D_n = \frac{(-1)^n H_0^{2n}}{4^n n! G_0^{2n}}. \qquad (70)$$

This expression was found by explicitly applying the recursive procedure - see Eqs.(63) and (65) - to the first terms and looking for a general formula, which must now be rigorously proved. To begin with, we note that the formula gives the correct value $D_0 = 1$, for $n = 0$. Now, we can calculate $D_{n+1}$ with the procedure of Section IIC to verify that the assumption made on $D_n$ is correct. To this aim, we firstly determine $\Lambda_{n+1}$ through Eq. (63) as follows

$$\Lambda_{n+1} = \frac{(-1)^{n+1} H_0^{2(n+1)}}{4^{n+1} n! G_0^{2(n+1)}}. \qquad (71)$$

Then, we can determine $D_{n+1}$ by means of Eq. (65)

$$\begin{align*}
D_{n+1} &= \int_y^x \frac{\mathcal{D}'(\eta)\mathcal{D}(\eta, y)^n}{\mathcal{D}(x, y)^{n+1}} \frac{(-1)^{n+1} H_0^{2(n+1)}}{4^{n+1}n!G_0^{2(n+1)}} d\eta \\
&= \int_y^x \frac{(\eta - y)^n}{(x - y)^{n+1}} \frac{(-1)^{n+1} H_0^{2(n+1)}}{4^{n+1}n!G_0^{2(n+1)}} d\eta \\
&= \frac{(-1)^{n+1} H_0^{2(n+1)}}{4^{n+1}(n+1)!G_0^{2(n+1)}}, \tag{72}
\end{align*}$$

which is consistent with Eq. (70) through the replacement $n \to n+1$. This confirms the validity of the result by mathematical induction.

This approach, based on an *ansatz* like that in Eq. (70), which is then verified by induction on the natural numbers is interesting but can probably be applied to a variety of special cases but not to the general problem. In fact, in the most complicated

cases it is difficult to propose an expression for $D_n$ by applying only a few recursions of the procedure. However, it remains useful for various applications discussed below.

We can now explicitly determine $F$ as follows

$$\begin{align*}
F(x,t;y,t_0) &= \sum_{n=0}^{+\infty} F_n(x;y)(t-t_0)^n \\
&= \sum_{n=0}^{+\infty} D_n(x)(t-t_0)^n F_0(x;y) \\
&= \sum_{n=0}^{+\infty} \frac{\left(-\frac{H_0^2(t-t_0)}{4G_0^2}\right)^n}{n!} F_0(x;y) \\
&= \exp\left[-\frac{H_0^2(t-t_0)}{4G_0^2}\right] \exp\left[-\frac{H_0}{2G_0^2}(x-y)\right]. \tag{73}
\end{align*}$$

We can finally multiply $K_0$ by $F$, and it is easily seen that we get the correct solution for $K$, given by Eq. (67). The series defining the function $F$ is clearly convergent and thus determines the first exponential function in Eq. (73). Although rather simple, this example shows that the coefficients $D_n$ decrease rapidly with $n$ and thus a few terms are often sufficient to have a good approximation on short-time intervals.

An illustration of this behaviour is found in Fig. 1, where we show $K_0$, $K$, and different approximations given by $K_0Q_m$, where $Q_m = \sum_{n=0}^m F_n(x;y)(t-t_0)^n$; $Q_m$ is partial sum of Eq. (73). We see that the first approximation $K_0Q_0$ is too crude and in fact is not even properly normalizable. On the other hand, the subsequent approximations are automatically normalized to within a negligible error. To conclude, we finally remark that the results of this first example are independent of $\alpha$ since the noise is additive and not multiplicative. We next turn to a first example with multiplicative noise.

## B. Geometric Brownian processes

Geometric Brownian processes find application in physics, e.g., in the statistical theory of turbulence [38, 39], but they are much more popular in finance, in particular in the modeling of stock option prices [53–56].

For this example we will see that we can compute the kernel with our recursive procedure for an arbitrary interpretation of the stochastic integrals (i.e. for any value of $\alpha$). We consider an exemplary geometric Brownian process described by the stochastic differential equation

$$\frac{dx}{dt} = -H_0x + G_0|x|\xi(t), \qquad (74)$$

where we considered the absolute value in the diffusion term in order to have $g(x) \ge 0$, $x \in \mathbb{R}$, that is a key assumption of all our work (considering $G_0 > 0$). We have $h(x) = -H_0x$, $g(x) = G_0|x|$, and $g'(x) = G_0 \operatorname{sgn}(x)$, where $\operatorname{sgn}(\zeta)$ represents the signum function returning the sign of $\zeta$, i.e., $\operatorname{sgn}(\zeta) = +1$ if $\zeta > 0$, $\operatorname{sgn}(\zeta) = 0$ if $\zeta = 0$, and $\operatorname{sgn}(\zeta) = -1$ if $\zeta < 0$. It
<figure><img src="image_2.png" /><figcaption>FIGURE 2. Convergence of the short-time expansion for a geometric Brownian process. We plot the short-time propagator $K_0$, the full propagator $K$, and some low-order approximations $K_0Q_0$, $K_0Q_1$, and $K_0Q_2$, to illustrate the convergence $K_0Q_m \to K$ when $m \to \infty$. We adopted the parameter values $H_0 = -1$, $G_0 = 1$, $y = 5$, $t_0 = 0$, $t = 2$, and $\alpha = 1/2$ (in arbitrary units).</figcaption></figure>

easily follows from Eqs.(5), (6), and (58) that

$$f(x) = -H_0 + 2(\alpha - 1)G_0^2, \quad (75)$$

$$\ell(x) = [-H_0 + 2(\alpha - 2)G_0^2] x, \quad (76)$$

$$\Phi(x) = \left[-\frac{H_0}{G_0^2} + (2\alpha - 3)\right] \frac{1}{2x}, \quad (77)$$

valid for any $x \in \mathbb{R}$. An important point concerns the determination of the function $\mathcal{D}(x,y)$, defined as

$$\mathcal{D}(x,y) = \int_y^x \frac{d\eta}{g(\eta)} = \frac{1}{G_0} \int_y^x \frac{d\eta}{|\eta|}. \quad (78)$$

We know that $\int d\eta/|\eta| = \operatorname{sgn}(\eta) \ln|\eta| + C$ (for some real constant $C$), but we must be careful in the use of this indefinite integral. In fact, if $x$ and $y$ have the same sign (both positive or both negative), there is no problem and the result can be written without doubt as $\int_y^x d\eta/|\eta| = \operatorname{sgn}(x) \ln|x/y| = \operatorname{sgn}(y) \ln|x/y|$. However, when $\operatorname{sgn}(x) \neq \operatorname{sgn}(y)$, the integration region certainly crosses the origin and the definite integral is then divergent, i.e. $\int_y^x d\eta/|\eta| \to \pm\infty$. To summarize, we have:

$$\mathcal{D}(x,y) = \begin{cases} \frac{1}{G_0} \operatorname{sgn}(x) \ln \left| \frac{x}{y} \right| & \text{if } \operatorname{sgn}(x) = \operatorname{sgn}(y), \\ \pm\infty & \text{if } \operatorname{sgn}(x) \neq \operatorname{sgn}(y). \end{cases} \quad (79)$$

This complex behavior of the function $\mathcal{D}$ has important implications for the short-time propagator structure. To begin, we can write

$$K_0(x,t;y,t_0) = \frac{1}{G_0|y|\sqrt{4\pi(t-t_0)}} \exp\left[-\frac{\mathcal{D}^2(x,y)}{4(t-t_0)}\right], \quad (80)$$

and therefore, we obtain

$$K_0(x,t;y,t_0) = \frac{\exp\left[-\frac{\ln^2\left|\frac{x}{y}\right|}{4G_0^2(t-t_0)}\right]}{G_0|y|\sqrt{4\pi(t-t_0)}} \mathbf{1}(x), \quad y > 0, \quad (81)$$

$$K_0(x,t;y,t_0) = \frac{\exp\left[-\frac{\ln^2\left|\frac{x}{y}\right|}{4G_0^2(t-t_0)}\right]}{G_0|y|\sqrt{4\pi(t-t_0)}} \mathbf{1}(-x), \quad y < 0, \quad (82)$$

where the Heaviside step function $\mathbf{1}(\zeta)$ is defined as $\mathbf{1}(\zeta) = 1$ if $\zeta \ge 0$, and $\mathbf{1}(\zeta) = 0$ if $\zeta < 0$. This means that the propagator vanishes on the negative real semi-axis when the initial value $y$ is positive, and vice versa, it vanishes on the positive real semi-axis, when the initial value $y$ is negative. This result is derived from the fact that the function $\mathcal{D}$ is infinite when the signs of $x$ and $y$ are discordant, and then we observe that $K_0 \to 0$ in this case. Therefore, the stochastic process cannot cross the point $x = 0$, as if there were a fictitious reflecting boundary condition. This is consistent with the fact that both the drift and diffusion terms are zero for $x = 0$.

The short-time propagator must be now improved by means of the correction function $F$. As before, we apply the recursive procedure stated in Eqs.(63) and (65) and the observation of the first terms leads to the expression

$$D_n = \frac{(-1)^n G_0^{2n}}{4^n n!} \left(2\alpha - 1 - \frac{H_0}{G_0^2}\right)^{2n}, \quad (83)$$

which must be proved by mathematical induction. Firstly, for $n = 0$, we correctly have $D_0 = 1$. Based on Eq. (83), the straightforward application of Eq. (63) delivers the following coefficient $\Lambda_{n+1}$

$$\Lambda_{n+1} = \frac{(-1)^{n+1} G_0^{2(n+1)}}{4^{n+1} n!} \left(2\alpha - 1 - \frac{H_0}{G_0^2}\right)^{2(n+1)}. \quad (84)$$

For applying Eq. (65), we observe that $\Lambda_{n+1}$ is a constant, and we can restrict the study the case with $\operatorname{sgn}(x) = \operatorname{sgn}(y)$ for what has already been discussed in the calculation of $K_0$ (indeed, $K_0$ is zero if $\operatorname{sgn}(x) \neq \operatorname{sgn}(y)$). We then write

$$\begin{aligned} D_{n+1} &= \frac{\Lambda_{n+1}}{[\operatorname{sgn}(x) \ln |\frac{x}{y}|]^{n+1}} \int_y^x \frac{[\operatorname{sgn}(\eta) \ln |\frac{\eta}{y}|]^n}{|\eta|} d\eta \\ &= \frac{\Lambda_{n+1}}{n+1} \\ &= \frac{(-1)^{n+1} G_0^{2(n+1)}}{4^{n+1}(n+1)!} \left(2\alpha - 1 - \frac{H_0}{G_0^2}\right)^{2(n+1)}, \quad (85) \end{aligned}$$

which is consistent with Eq. (83) as a result of the substitution $n \to n+1$. The mathematical induction finally proves the expression for $D_n$, which can be used to obtain the correction
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Convergence of the short-time expansion for a geometric Brownian process. We plot the short-time propagator $K_0$, the full propagator $K$, and some low-order approximations $K_0Q_0$, $K_0Q_1$, and $K_0Q_2$, to illustrate the convergence $K_0Q_m \to K$ when $m \to \infty$. We adopted the parameter values $H_0 = -1$, $G_0 = 1$, $y = 5$, $t_0 = 0$, $t = 2$, and $\alpha = 1/2$ (in arbitrary units).</figcaption></figure>

easily follows from Eqs.(5), (6), and (58) that

$$f(x) = -H_0 + 2(\alpha - 1)G_0^2, \quad (75)$$

$$\ell(x) = [-H_0 + 2(\alpha - 2)G_0^2] x, \quad (76)$$

$$\Phi(x) = \left[ -\frac{H_0}{G_0^2} + (2\alpha - 3) \right] \frac{1}{2x}, \quad (77)$$

valid for any $x \in \mathbb{R}$. An important point concerns the determination of the function $\mathcal{D}(x,y)$, defined as

$$\mathcal{D}(x,y) = \int_y^x \frac{d\eta}{g(\eta)} = \frac{1}{G_0} \int_y^x \frac{d\eta}{|\eta|}. \quad (78)$$

We know that $\int d\eta/|\eta| = \operatorname{sgn}(\eta) \ln|\eta| + C$ (for some real constant $C$), but we must be careful in the use of this indefinite integral. In fact, if $x$ and $y$ have the same sign (both positive or both negative), there is no problem and the result can be written without doubt as $\int_y^x d\eta/|\eta| = \operatorname{sgn}(x) \ln|x/y| = \operatorname{sgn}(y) \ln|x/y|$. However, when $\operatorname{sgn}(x) \neq \operatorname{sgn}(y)$, the integration region certainly crosses the origin and the definite integral is then divergent, i.e. $\int_y^x d\eta/|\eta| \to \pm\infty$. To summarize, we have:

$$\mathcal{D}(x,y) = \begin{cases} \frac{1}{G_0} \operatorname{sgn}(x) \ln \left| \frac{x}{y} \right| & \text{if } \operatorname{sgn}(x) = \operatorname{sgn}(y), \\ \pm\infty & \text{if } \operatorname{sgn}(x) \neq \operatorname{sgn}(y). \end{cases} \quad (79)$$

This complex behavior of the function $\mathcal{D}$ has important implications for the short-time propagator structure. To begin, we can write

$$K_0(x,t;y,t_0) = \frac{1}{G_0|y|\sqrt{4\pi(t-t_0)}} \exp\left[-\frac{\mathcal{D}^2(x,y)}{4(t-t_0)}\right], \quad (80)$$

and therefore, we obtain

$$K_0(x,t;y,t_0) = \frac{\exp\left[-\frac{\ln^2\left|\frac{x}{y}\right|}{4G_0^2(t-t_0)}\right]}{G_0|y|\sqrt{4\pi(t-t_0)}} \mathbf{1}(x), \quad y > 0, \quad (81)$$

$$K_0(x,t;y,t_0) = \frac{\exp\left[-\frac{\ln^2\left|\frac{x}{y}\right|}{4G_0^2(t-t_0)}\right]}{G_0|y|\sqrt{4\pi(t-t_0)}} \mathbf{1}(-x), \quad y < 0, \quad (82)$$

where the Heaviside step function $\mathbf{1}(\zeta)$ is defined as $\mathbf{1}(\zeta) = 1$ if $\zeta \ge 0$, and $\mathbf{1}(\zeta) = 0$ if $\zeta < 0$. This means that the propagator vanishes on the negative real semi-axis when the initial value $y$ is positive, and vice versa, it vanishes on the positive real semi-axis, when the initial value $y$ is negative. This result is derived from the fact that the function $\mathcal{D}$ is infinite when the signs of $x$ and $y$ are discordant, and then we observe that $K_0 \to 0$ in this case. Therefore, the stochastic process cannot cross the point $x = 0$, as if there were a fictitious reflecting boundary condition. This is consistent with the fact that both the drift and diffusion terms are zero for $x = 0$.

The short-time propagator must be now improved by means of the correction function $F$. As before, we apply the recursive procedure stated in Eqs.(63) and (65) and the observation of the first terms leads to the expression

$$D_n = \frac{(-1)^n G_0^{2n}}{4^n n!} \left( 2\alpha - 1 - \frac{H_0}{G_0^2} \right)^{2n}, \quad (83)$$

which must be proved by mathematical induction. Firstly, for $n = 0$, we correctly have $D_0 = 1$. Based on Eq. (83), the straightforward application of Eq. (63) delivers the following coefficient $\Lambda_{n+1}$

$$\Lambda_{n+1} = \frac{(-1)^{n+1} G_0^{2(n+1)}}{4^{n+1} n!} \left( 2\alpha - 1 - \frac{H_0}{G_0^2} \right)^{2(n+1)}. \quad (84)$$

For applying Eq. (65), we observe that $\Lambda_{n+1}$ is a constant, and we can restrict the study the case with $\operatorname{sgn}(x) = \operatorname{sgn}(y)$ for what has already been discussed in the calculation of $K_0$ (indeed, $K_0$ is zero if $\operatorname{sgn}(x) \neq \operatorname{sgn}(y)$). We then write

$$\begin{aligned} D_{n+1} &= \frac{\Lambda_{n+1}}{[\operatorname{sgn}(x) \ln |\frac{x}{y}|]^{n+1}} \int_y^x \frac{[\operatorname{sgn}(\eta) \ln |\frac{\eta}{y}|]^n}{|\eta|} d\eta \\ &= \frac{\Lambda_{n+1}}{n+1} \\ &= \frac{(-1)^{n+1} G_0^{2(n+1)}}{4^{n+1}(n+1)!} \left( 2\alpha - 1 - \frac{H_0}{G_0^2} \right)^{2(n+1)}, \quad (85) \end{aligned}$$

which is consistent with Eq. (83) as a result of the substitution $n \to n+1$. The mathematical induction finally proves the expression for $D_n$, which can be used to obtain the correction
$F = F_0 \sum_{n=0}^{+\infty} D_n(t-t_0)^n$. We determine

$$F_0 = \exp \left[ \int_y^x \Phi(\eta) d\eta \right] = \exp \left[ \left( -\frac{H_0}{2G_0^2} + \frac{2\alpha - 3}{2} \right) \ln \left| \frac{x}{y} \right| \right], \quad (86)$$

and

$$\begin{aligned} F &= F_0 \sum_{n=0}^{+\infty} \frac{(-1)^n G_0^{2n}}{4^n n!} \left( 2\alpha - 1 - \frac{H_0}{G_0^2} \right)^{2n} (t-t_0)^n \\ &= F_0 \exp \left[ -\frac{G_0^2}{4} \left( 2\alpha - 1 - \frac{H_0}{G_0^2} \right)^2 (t-t_0) \right]. \end{aligned} \quad (87)$$

To conclude, after straightforward calculations we obtain the complete propagator $K(x, t; y, t_0) = K_0 F$ in the form

$$K = \frac{\exp \left[ - \frac{\left[ \ln \left| \frac{x}{y} \right| + H_0(t-t_0) - (2\alpha-1)G_0^2(t-t_0) \right]^2}{4G_0^2(t-t_0)} \right]}{G_0 |x| \sqrt{4\pi(t-t_0)}} \mathbf{1}(x), \quad y > 0, \quad (88)$$

$$K = \frac{\exp \left[ - \frac{\left[ \ln \left| \frac{x}{y} \right| + H_0(t-t_0) - (2\alpha-1)G_0^2(t-t_0) \right]^2}{4G_0^2(t-t_0)} \right]}{G_0 |x| \sqrt{4\pi(t-t_0)}} \mathbf{1}(-x), \quad y < 0, \quad (89)$$

which is in full agreement with previous investigations [16, 57]. A numerical result is displayed in Fig. 2, where we show the fast convergence of the process. We plot the short-time propagator $K_0$, the full propagator $K$, and some approximations $K_0 Q_0$, $K_0 Q_1$, and $K_0 Q_2$, to show the convergence $K_0 Q_m \to K$, when $m \to \infty$. Here, $Q_m = \sum_{n=0}^m F_n(x; y)(t-t_0)^n$, representing the partial sum of Eq. (87). We have rigorously proved by our method that if the geometric stochastic equation is defined on the entire real axis, then the support of the propagator consists of only the semi-axis in which the initial value is present. The obtained result for the propagator is the log-normal distribution which is typically associated with geometric Brownian motion. We note that the result retains the dependence on the discretization parameter $\alpha$ and thus holds for all definitions of the the stochastic integral. However, it may be worth noting that in most applications the geometric Brownian process is considered in the Itô-interpretation, with $\alpha = 0$.

### C. A stochastic process for chromatin remodeling

We will now introduce two practical applications for which the full propagator is difficult to obtain, therefore the approximation of the short-term expression of $K$ may have an important, heuristic value.

We firstly consider a stochastic model that has been used to map the stochastic motion of a chromatin remodeler acting on a nucleosome to an active Brownian dimer [26–28]. Originally developed for general dimeric motors, it has recently been linked via structural biology insights on the ac-

tion of two motor domains of a chromatin remodeling enzyme on the DNA wrapped around a nucleosome. In this understanding, first one of the motor domain pulls in DNA into the nucleosome, and subsequently the second motor domain expels the surplus DNA outward, thereby establishing a proper wrapping of 147 base pairs of DNA around the protein core of the nucleosome. The motion of the two-domain nucleosome-remodeler complex can then be separated into a center of mass motion, and a stochastic differential equation for the relative motion, which is given by a very simple expression in view of the complexity of the system

$$\frac{dx}{dt} = h(x) + g(x)\xi(t) \quad (90)$$

with the drift and diffusion terms given by

$$h(x) = - \left( \frac{1}{\beta} + \frac{1}{\beta + \gamma x} \right) kx, \quad (91)$$

$$g(x) = \frac{A}{(\beta + \gamma x)^2}. \quad (92)$$

These expressions are valid when only active, ATP-dependent fluctuations are retained and thermal fluctuations are ignored. The equation is studied under the Fisk-Stratonovich interpretation ($\alpha = 1/2$). This choice was made only to fix ideas, and any other value for $\alpha$ could have been considered. In physical and biophysical applications, however, this is the value typically adopted.

The parameter $A$ signifies the strength of the active process, while the constants $\beta$ and $\gamma$ are friction terms ($\gamma \ll \beta$); $k$ is the spring constant of a harmonic interaction potential between the remodeler motor domains on DNA. We stress that this model is purely phenomenological and therefore also other forms of these functions can be adopted and can in fact be found in the literature [26, 27].

Applying our formal analysis, we fix the initial position $y = 0$, and the initial time $t_0 = 0$. First of all, we determine the quantity $\mathcal{D}$ as

$$\mathcal{D} = \frac{1}{A} \left( \frac{1}{3} \gamma^2 x^3 + \beta \gamma x^2 + \beta^2 x \right), \quad (93)$$

and therefore we obtain the short-time propagator in the form

$$K_0(x, t; y, t_0) = \beta^2 \frac{\exp \left[ - \frac{1}{4A^2t} \left( \frac{1}{3}\gamma^2 x^3 + \beta\gamma x^2 + \beta^2 x \right)^2 \right]}{\sqrt{4\pi A^2t}}. \quad (94)$$

Moreover, the first correction $F_0$ is found as

$$F_0 = \left( \frac{\gamma x + \beta}{\beta} \right)^2 \exp \left( - \frac{\beta^3 k x^2}{2A^2} - \frac{7\beta^2 \gamma k x^3}{6A^2} - \frac{9\beta \gamma^2 k x^4}{8A^2} - \frac{\gamma^3 k x^5}{2A^2} - \frac{\gamma^4 k x^6}{12\beta A^2} \right). \quad (95)$$

The application of the recursive procedure given by Eqs. (63)
$F = F_0 \sum_{n=0}^{+\infty} D_n(t-t_0)^n$. We determine

$$F_0 = \exp \left[ \int_y^x \Phi(\eta) d\eta \right] = \exp \left[ \left( -\frac{H_0}{2G_0^2} + \frac{2\alpha - 3}{2} \right) \ln \left| \frac{x}{y} \right| \right], \quad (86)$$

and

$$ \begin{aligned} F &= F_0 \sum_{n=0}^{+\infty} \frac{(-1)^n G_0^{2n}}{4^n n!} \left( 2\alpha - 1 - \frac{H_0}{G_0^2} \right)^{2n} (t-t_0)^n \\ &= F_0 \exp \left[ -\frac{G_0^2}{4} \left( 2\alpha - 1 - \frac{H_0}{G_0^2} \right)^2 (t-t_0) \right]. \end{aligned} \quad (87) $$

To conclude, after straightforward calculations we obtain the complete propagator $K(x, t; y, t_0) = K_0 F$ in the form

$$K = \frac{\exp \left[ - \frac{\left[ \ln \left| \frac{x}{y} \right| + H_0 (t-t_0) - (2\alpha-1)G_0^2 (t-t_0) \right]^2}{4G_0^2 (t-t_0)} \right]}{G_0 |x| \sqrt{4\pi(t-t_0)}} \mathbf{1}(x), \quad y > 0, \quad (88)$$

$$K = \frac{\exp \left[ - \frac{\left[ \ln \left| \frac{x}{y} \right| + H_0 (t-t_0) - (2\alpha-1)G_0^2 (t-t_0) \right]^2}{4G_0^2 (t-t_0)} \right]}{G_0 |x| \sqrt{4\pi(t-t_0)}} \mathbf{1}(-x), \quad y < 0, \quad (89)$$

which is in full agreement with previous investigations [16, 57]. A numerical result is displayed in Fig. 2, where we show the fast convergence of the process. We plot the short-time propagator $K_0$, the full propagator $K$, and some approximations $K_0 Q_0$, $K_0 Q_1$, and $K_0 Q_2$, to show the convergence $K_0 Q_m \to K$, when $m \to \infty$. Here, $Q_m = \sum_{n=0}^m F_n(x; y)(t-t_0)^n$, representing the partial sum of Eq. (87). We have rigorously proved by our method that if the geometric stochastic equation is defined on the entire real axis, then the support of the propagator consists of only the semi-axis in which the initial value is present. The obtained result for the propagator is the log-normal distribution which is typically associated with geometric Brownian motion. We note that the result retains the dependence on the discretization parameter $\alpha$ and thus holds for all definitions of the the stochastic integral. However, it may be worth noting that in most applications the geometric Brownian process is considered in the Itô-interpretation, with $\alpha = 0$.

### C. A stochastic process for chromatin remodeling

We will now introduce two practical applications for which the full propagator is difficult to obtain, therefore the approximation of the short-term expression of $K$ may have an important, heuristic value.

We firstly consider a stochastic model that has been used to map the stochastic motion of a chromatin remodeler acting on a nucleosome to an active Brownian dimer [26–28]. Originally developed for general dimeric motors, it has recently been linked via structural biology insights on the ac-

tion of two motor domains of a chromatin remodeling enzyme on the DNA wrapped around a nucleosome. In this understanding, first one of the motor domain pulls in DNA into the nucleosome, and subsequently the second motor domain expels the surplus DNA outward, thereby establishing a proper wrapping of 147 base pairs of DNA around the protein core of the nucleosome. The motion of the two-domain nucleosome-remodeler complex can then be separated into a center of mass motion, and a stochastic differential equation for the relative motion, which is given by a very simple expression in view of the complexity of the system

$$\frac{dx}{dt} = h(x) + g(x) \xi(t) \quad (90)$$

with the drift and diffusion terms given by

$$h(x) = - \left( \frac{1}{\beta} + \frac{1}{\beta + \gamma x} \right) kx, \quad (91)$$

$$g(x) = \frac{A}{(\beta + \gamma x)^2}. \quad (92)$$

These expressions are valid when only active, ATP-dependent fluctuations are retained and thermal fluctuations are ignored. The equation is studied under the Fisk-Stratonovich interpretation ($\alpha = 1/2$). This choice was made only to fix ideas, and any other value for $\alpha$ could have been considered. In physical and biophysical applications, however, this is the value typically adopted.

The parameter $A$ signifies the strength of the active process, while the constants $\beta$ and $\gamma$ are friction terms ($\gamma \ll \beta$); $k$ is the spring constant of a harmonic interaction potential between the remodeler motor domains on DNA. We stress that this model is purely phenomenological and therefore also other forms of these functions can be adopted and can in fact be found in the literature [26, 27].

Applying our formal analysis, we fix the initial position $y = 0$, and the initial time $t_0 = 0$. First of all, we determine the quantity $\mathcal{D}$ as

$$\mathcal{D} = \frac{1}{A} \left( \frac{1}{3} \gamma^2 x^3 + \beta \gamma x^2 + \beta^2 x \right), \quad (93)$$

and therefore we obtain the short-time propagator in the form

$$K_0(x, t; y, t_0) = \beta^2 \frac{\exp \left[ - \frac{1}{4A^2t} \left( \frac{1}{3}\gamma^2 x^3 + \beta\gamma x^2 + \beta^2 x \right)^2 \right]}{\sqrt{4\pi A^2t}}. \quad (94)$$

Moreover, the first correction $F_0$ is found as

$$F_0 = \left( \frac{\gamma x + \beta}{\beta} \right)^2 \exp \left( - \frac{\beta^3 k x^2}{2A^2} - \frac{7\beta^2 \gamma k x^3}{6A^2} - \frac{9\beta\gamma^2 k x^4}{8A^2} - \frac{\gamma^3 k x^5}{2A^2} - \frac{\gamma^4 k x^6}{12\beta A^2} \right). \quad (95)$$

The application of the recursive procedure given by Eqs.(63)
<figure><img src="image_3.png" /><figcaption>FIGURE 3. Application of the recursive procedure to the stochastic equation for the remodeler-nucleosome complex. In the left panel we adopted the parameters $k = 2$, $\beta = 1$, $\gamma = 1/10$, and $A = 1$, and in the right panel $k = 2$, $\beta = 1$, $\gamma = 1/10$, and $A = 2$ (in arbitrary units). In both cases we plot the short-time propagator $K_0$, and low-order approximations $K_0Q_0$, $K_0Q_1$, $K_0Q_2$, $K_0Q_3$, and $K_0Q_4$, to show the convergence of $K_0Q_m$ for $m \to \infty$. Here, $Q_m = \sum_{n=0}^{m} F_n(x;y)(t-t_0)^n$ represents the partial sum of Eq. (31). In both panels we adopted $t_0 = 0$, and $t = 1/2$.</figcaption></figure>

and (65) yields closed form expressions for the coefficients $D_n$. While these quite lengthy expressions can be easily obtained through standard symbolic environments like MAPLE, we display in the Appendix only the series expansion of the coefficients $D_n$ with respect to the variable $x$, keeping fixed the other physical parameters.

The plot of the successive approximations for the propagator can be found in the left panel of Fig. 3 for a first set of parameters. Here we represented the functions $K_0$, and some approximations $K_0Q_0$, $K_0Q_1$, $K_0Q_2$, $K_0Q_3$, and $K_0Q_4$, to show the convergence of $K_0Q_m$ for $m \to \infty$. Here, $Q_m = \sum_{n=0}^{m} F_n(x;y)(t-t_0)^n$ represents the partial sum of Eq. (31). It means that $Q_0 = K_0$, $Q_1 = F_0(1+D_1t)$, $Q_2 = F_0(1+D_1t+D_2t^2)$, $Q_3 = F_0(1+D_1t+D_2t^2+D_3t^3)$, and $Q_4 = F_0(1+D_1t+D_2t^2+D_3t^3+D_4t^4)$. Here, we adopted $t_0 = 0$, and $t = 1/2$. It is interesting to note that a rather small number of terms is sufficient to have reasonably convergent behaviour. This can be seen from the fact that the last two approximations are almost overlapping and also from the fact that they are automatically well normalized. We emphasize that the first approximation $K_0F_0$ is again so poorly approximate that it is not even normalized.

The results for a second set of parameters, where we adopted a larger intensity of the noise related to the ATP-consumption, can be found in the right panel of Fig. 3. In this case, we observe a larger variance of the variable $x$, which induces larger propagator spread. Also in this case we can remark a good convergence with a quite limited number of terms. This means that the recursive procedure proposed in this paper can be effectively applied to physical and biological problems of some interest with relative simplicity. In particular, it makes it possible to obtain good approximations of the propagator of one-dimensional Fokker-Planck equations even when such equations are very difficult to solve by other approaches.

### D. Exponential heterogeneous diffusion

Our final example in this Section deals with the heterogeneous exponential diffusion described by $g(x) = G_0 \exp(\gamma x)$, with $G_0 > 0$ and $\gamma$ real. Such an exponential behavior of the diffusion coefficient has been used in a different biological context: it had been introduced to give a stochastic interpretation of the movement of the parasitic nematode *Phasmarhabditis hermaphrodita* [58–60]. Another biological example is given by an exponential rate of morphogen degradation in a reaction-subdiffusion model for cell development [61]. Morphogens are special signaling molecules whose spatial distribution controls the growth of embryonic cells [62]. However, such biologically-motivated examples are not the only ones for which an exponential diffusion coefficient can be used. For example, they were similarly employed to model the diffusion of impurities induced by irradiation [63, 64] : in this case, the exponential dependence originates from the decay of the radiation propagating in the material under study. Finally, they were also applied to the description of grain-boundary diffusion in nanostructured materials [65].

We thus consider the differential equation with exponential diffusion:

$$ \frac{dx}{dt} = G_0 \exp(\gamma x) \xi(t), \qquad (96) $$

where at first an arbitrary stochastic interpretation is assumed. We start our analysis by applying Eqs.(5), (6), and (58), which lead to the expressions

$$ f(x) = 4(\alpha - 1)G_0^2 \gamma^2 \exp(2\gamma x), \qquad (97) $$

$$ \ell(x) = 2(\alpha - 2)G_0^2 \gamma \exp(2\gamma x), \qquad (98) $$

$$ \Phi(x) = \frac{\gamma}{2}(2\alpha - 3). \qquad (99) $$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Application of the recursive procedure to the stochastic equation for the remodeler-nucleosome complex. In the left panel we adopted the parameters $k = 2$, $\beta = 1$, $\gamma = 1/10$, and $A = 1$, and in the right panel $k = 2$, $\beta = 1$, $\gamma = 1/10$, and $A = 2$ (in arbitrary units). In both cases we plot the short-time propagator $K_0$, and low-order approximations $K_0Q_0$, $K_0Q_1$, $K_0Q_2$, $K_0Q_3$, and $K_0Q_4$, to show the convergence of $K_0Q_m$ for $m \to \infty$. Here, $Q_m = \sum_{n=0}^{m} F_n(x;y)(t-t_0)^n$ represents the partial sum of Eq. (31). In both panels we adopted $t_0 = 0$, and $t = 1/2$.</figcaption></figure>

and (65) yields closed form expressions for the coefficients $D_n$. While these quite lengthy expressions can be easily obtained through standard symbolic environments like MAPLE, we display in the Appendix only the series expansion of the coefficients $D_n$ with respect to the variable $x$, keeping fixed the other physical parameters.

The plot of the successive approximations for the propagator can be found in the left panel of Fig. 3 for a first set of parameters. Here we represented the functions $K_0$, and some approximations $K_0Q_0$, $K_0Q_1$, $K_0Q_2$, $K_0Q_3$, and $K_0Q_4$, to show the convergence of $K_0Q_m$ for $m \to \infty$. Here, $Q_m = \sum_{n=0}^{m} F_n(x;y)(t-t_0)^n$ represents the partial sum of Eq. (31). It means that $Q_0 = K_0$, $Q_1 = F_0(1+D_1t)$, $Q_2 = F_0(1+D_1t+D_2t^2)$, $Q_3 = F_0(1+D_1t+D_2t^2+D_3t^3)$, and $Q_4 = F_0(1+D_1t+D_2t^2+D_3t^3+D_4t^4)$. Here, we adopted $t_0 = 0$, and $t = 1/2$. It is interesting to note that a rather small number of terms is sufficient to have reasonably convergent behaviour. This can be seen from the fact that the last two approximations are almost overlapping and also from the fact that they are automatically well normalized. We emphasize that the first approximation $K_0F_0$ is again so poorly approximate that it is not even normalized.

The results for a second set of parameters, where we adopted a larger intensity of the noise related to the ATP-consumption, can be found in the right panel of Fig. 3. In this case, we observe a larger variance of the variable $x$, which induces larger propagator spread. Also in this case we can remark a good convergence with a quite limited number of terms. This means that the recursive procedure proposed in this paper can be effectively applied to physical and biological problems of some interest with relative simplicity. In particular, it makes it possible to obtain good approximations of the propagator of one-dimensional Fokker-Planck equations even when such equations are very difficult to solve by other approaches.

### D. Exponential heterogeneous diffusion

Our final example in this Section deals with the heterogeneous exponential diffusion described by $g(x) = G_0 \exp(\gamma x)$, with $G_0 > 0$ and $\gamma$ real. Such an exponential behavior of the diffusion coefficient has been used in a different biological context: it had been introduced to give a stochastic interpretation of the movement of the parasitic nematode *Phasmarhabditis hermaphrodita* [58–60]. Another biological example is given by an exponential rate of morphogen degradation in a reaction-subdiffusion model for cell development [61]. Morphogens are special signaling molecules whose spatial distribution controls the growth of embryonic cells [62]. However, such biologically-motivated examples are not the only ones for which an exponential diffusion coefficient can be used. For example, they were similarly employed to model the diffusion of impurities induced by irradiation [63, 64] : in this case, the exponential dependence originates from the decay of the radiation propagating in the material under study. Finally, they were also applied to the description of grain-boundary diffusion in nanostructured materials [65].

We thus consider the differential equation with exponential diffusion:

$$\frac{dx}{dt} = G_0 \exp(\gamma x) \xi(t), \qquad (96)$$

where at first an arbitrary stochastic interpretation is assumed. We start our analysis by applying Eqs.(5), (6), and (58), which lead to the expressions

$$f(x) = 4(\alpha - 1)G_0^2 \gamma^2 \exp(2\gamma x), \qquad (97)$$

$$\ell(x) = 2(\alpha - 2)G_0^2 \gamma \exp(2\gamma x), \qquad (98)$$

$$\Phi(x) = \frac{\gamma}{2}(2\alpha - 3). \qquad (99)$$
As done previously, we applied the recursion for the first terms manually, see Eqs.(63) and (65), and realized that the parameters $D_n$ are described by the relation

$$D_n = \frac{(-1)^n G_0^{2n} \gamma^{2n} \exp[n\gamma(x+y)]}{4^n n!} \prod_{k=0}^{n-1} [4\alpha^2 - (2k+1)^2], \quad (100)$$

which we now demonstrate by mathematical induction. Firstly, for $n = 0$ we have $D_0 = 1$, considering that the empty product is unitary. Next, let's assume that the formula is true for a natural number $n$, and since $dD_n/dx = n\gamma D_n$, and $d\Phi/dx = 0$, we get from Eq. (63)

$$\Lambda_{n+1} = -G_0^2 \gamma^2 \exp(2\gamma x) D_n \frac{4\alpha^2 - (2n+1)^2}{4}. \quad (101)$$

We introduce the polynomials $Q_n(\zeta)$ such that $Q_0(\zeta) = 1$ and $Q_{n+1}(\zeta) = [4\zeta^2 - (2n+1)^2]Q_n(\zeta)$. They allow us to write Eq. (100) in the form

$$D_n = \frac{(-1)^n G_0^{2n} \gamma^{2n} \exp[n\gamma(x+y)]}{4^n n!} Q_n(\alpha), \quad (102)$$

Now, by means of the expression $\mathcal{D} = [\exp(-\gamma y) - \exp(-\gamma x)]/G_0\gamma$, we can apply Eq. (65) and we obtain after straightforward calculation

$$D_{n+1} = \frac{(-1)^{n+1} G_0^{2(n+1)} \gamma^{2(n+1)} \exp[(n+1)\gamma(x+y)]}{4^{n+1}(n+1)!} Q_{n+1}(\alpha), \quad (103)$$

which finally proves Eq. (100) by mathematical induction.

In this case we have obtained the closed form of the Taylor coefficients but it is important now to remark that the corresponding series leads to a meaningful result only if $\alpha = 1/2$. In fact, with $\alpha = 1/2$, all coefficients $D_n$ vanish with $n \ge 1$ as seen by observing the product in Eq. (100). The expansion is thus cut off after the first term. On the other hand, it is seen that with $\alpha \ne 1/2$ the series of powers with coefficients $D_n$ cannot converge because of the divergent behavior of $D_n$ with $n$ since the product grows roughly as $(n!)^2$. Therefore, only the Fisk-Stratonovich interpretation leads to a well-defined development for the exponential diffusion (this case will be further discussed in the next Section). As a consequence the probability density for the exponential diffusive process cannot be developed in a Taylor series independently from the interpretation of the stochastic integrals.

Although it is interesting to have understood that the solution by series can be obtained only for certain values of $\alpha$, a theoretical understanding of this fact is lacking. In particular, for the time being it is impossible to predict for a certain problem what are the values of $\alpha$ for which we have a convergent solution. This point certainly requires further theoretical analysis in the near future.

## IV. STOCHASTIC PROCESSES IN THE FISK-STRATONOVICH INTERPRETATION - NEW INSIGHTS FROM OUR FORMALISM

In this Section we take a different point of view. Rather than starting out from a specifically defined model, as we did in Section III, we here develop rather general consequences following from of our formalism. In this Section we consider stochastic processes always in the Fisk-Stratonovich calculus. We first report on some general results in the case of pure diffusions, i.e. in the absence of a drift term. Then we turn to the special case in which drift and diffusion terms conspire in such a way as to render our recursion equations particularly simple to solve.

### A. The case of pure diffusion and no drift

We first consider a stochastic process without drift, $h(x) = 0$, interpreted through the Fisk-Stratonovich integration with $\alpha = 1/2$. We have

$$\frac{dx}{dt} = g(x)\xi(t), \quad (104)$$

with $g(x) \ge 0$ for any $x \in \mathbb{R}$. To take advantage of our procedure, we can directly apply Eqs.(5), (6), and (58) and obtain the expressions

$$f(x) = -g'^2(x) - g(x)g''(x), \quad (105)$$

$$\ell(x) = -3g(x)g'(x), \quad (106)$$

$$\Phi(x) = -\frac{g'(x)}{g(x)}. \quad (107)$$

In order to implement the recursion, we start form $D_0 = 1$ and we calculate $\Lambda_1$ through Eq. (63), eventually obtaining

$$\begin{aligned} \Lambda_1(x) &= g'^2(x) + g(x)g''(x) - 3g(x)g'(x)\frac{g'(x)}{g(x)} \\ &\quad + g^2(x)\left(-\frac{g''(x)g(x) - g'^2(x)}{g^2(x)} + \frac{g'^2(x)}{g^2(x)}\right) \\ &= g'^2(x) + g(x)g''(x) - 3g'^2(x) - g''(x)g(x) \\ &\quad + g'^2(x) + g^2(x) = 0. \end{aligned} \quad (108)$$

The fact that $\Lambda_1 = 0$, by induction, implies that all parameters $\Lambda_n$ and $D_n$ are zero from $n = 1$ until infinity. Hence, the correction function $F$ is simply given by the first term $F_0$, which can be calculated as follows

$$F_0(x; y) = \exp \left[ \int_y^x \Phi(\eta, y) d\eta \right] = \frac{g(y)}{g(x)}, \quad (109)$$

where we used Eq. (107). To conclude, we obtain the closed form expression of the propagator for Eq. (104) as $K = K_0 F_0$,
As done previously, we applied the recursion for the first terms manually, see Eqs.(63) and (65), and realized that the parameters $D_n$ are described by the relation

$$D_n = \frac{(-1)^n G_0^{2n} \gamma^{2n} \exp[n\gamma(x+y)]}{4^n n!} \prod_{k=0}^{n-1} [4\alpha^2 - (2k+1)^2], \quad (100)$$

which we now demonstrate by mathematical induction. Firstly, for $n = 0$ we have $D_0 = 1$, considering that the empty product is unitary. Next, let's assume that the formula is true for a natural number $n$, and since $dD_n/dx = n\gamma D_n$, and $d\Phi/dx = 0$, we get from Eq. (63)

$$\Lambda_{n+1} = -G_0^2 \gamma^2 \exp(2\gamma x) D_n \frac{4\alpha^2 - (2n+1)^2}{4}. \quad (101)$$

We introduce the polynomials $Q_n(\zeta)$ such that $Q_0(\zeta) = 1$ and $Q_{n+1}(\zeta) = [4\zeta^2 - (2n+1)^2]Q_n(\zeta)$. They allow us to write Eq. (100) in the form

$$D_n = \frac{(-1)^n G_0^{2n} \gamma^{2n} \exp[n\gamma(x+y)]}{4^n n!} Q_n(\alpha), \quad (102)$$

Now, by means of the expression $\mathcal{D} = [\exp(-\gamma y) - \exp(-\gamma x)]/G_0\gamma$, we can apply Eq. (65) and we obtain after straightforward calculation

$$D_{n+1} = \frac{(-1)^{n+1} G_0^{2(n+1)} \gamma^{2(n+1)} \exp[(n+1)\gamma(x+y)]}{4^{n+1}(n+1)!} Q_{n+1}(\alpha), \quad (103)$$

which finally proves Eq. (100) by mathematical induction.

In this case we have obtained the closed form of the Taylor coefficients but it is important now to remark that the corresponding series leads to a meaningful result only if $\alpha = 1/2$. In fact, with $\alpha = 1/2$, all coefficients $D_n$ vanish with $n \ge 1$ as seen by observing the product in Eq. (100). The expansion is thus cut off after the first term. On the other hand, it is seen that with $\alpha \ne 1/2$ the series of powers with coefficients $D_n$ cannot converge because of the divergent behavior of $D_n$ with $n$ since the product grows roughly as $(n!)^2$. Therefore, only the Fisk-Stratonovich interpretation leads to a well-defined development for the exponential diffusion (this case will be further discussed in the next Section). As a consequence the probability density for the exponential diffusive process cannot be developed in a Taylor series independently from the interpretation of the stochastic integrals.

Although it is interesting to have understood that the solution by series can be obtained only for certain values of $\alpha$, a theoretical understanding of this fact is lacking. In particular, for the time being it is impossible to predict for a certain problem what are the values of $\alpha$ for which we have a convergent solution. This point certainly requires further theoretical analysis in the near future.

## IV. STOCHASTIC PROCESSES IN THE FISK-STRATONOVICH INTERPRETATION - NEW INSIGHTS FROM OUR FORMALISM

In this Section we take a different point of view. Rather than starting out from a specifically defined model, as we did in Section III, we here develop rather general consequences following from of our formalism. In this Section we consider stochastic processes always in the Fisk-Stratonovich calculus. We first report on some general results in the case of pure diffusions, i.e. in the absence of a drift term. Then we turn to the special case in which drift and diffusion terms conspire in such a way as to render our recursion equations particularly simple to solve.

### A. The case of pure diffusion and no drift

We first consider a stochastic process without drift, $h(x) = 0$, interpreted through the Fisk-Stratonovich integration with $\alpha = 1/2$. We have

$$\frac{dx}{dt} = g(x)\xi(t), \quad (104)$$

with $g(x) \ge 0$ for any $x \in \mathbb{R}$. To take advantage of our procedure, we can directly apply Eqs.(5), (6), and (58) and obtain the expressions

$$f(x) = -g'^2(x) - g(x)g''(x), \quad (105)$$

$$\ell(x) = -3g(x)g'(x), \quad (106)$$

$$\Phi(x) = -\frac{g'(x)}{g(x)}. \quad (107)$$

In order to implement the recursion, we start form $D_0 = 1$ and we calculate $\Lambda_1$ through Eq. (63), eventually obtaining

$$\begin{aligned} \Lambda_1(x) &= g'^2(x) + g(x)g''(x) - 3g(x)g'(x)\frac{g'(x)}{g(x)} \\ &\quad +g^2(x)\left(-\frac{g''(x)g(x)-g'^2(x)}{g^2(x)} + \frac{g'^2(x)}{g^2(x)}\right) \\ &= g'^2(x) + g(x)g''(x) - 3g'^2(x) - g''(x)g(x) \\ &\quad +g'^2(x) + g'^2(x) = 0. \end{aligned} \quad (108)$$

The fact that $\Lambda_1 = 0$, by induction, implies that all parameters $\Lambda_n$ and $D_n$ are zero from $n = 1$ until infinity. Hence, the correction function $F$ is simply given by the first term $F_0$, which can be calculated as follows

$$F_0(x;y) = \exp\left[\int_y^x \Phi(\eta,y)d\eta\right] = \frac{g(y)}{g(x)}, \quad (109)$$

where we used Eq. (107). To conclude, we obtain the closed form expression of the propagator for Eq. (104) as $K = K_0F_0$,
<figure><img src="image_4.png" /><figcaption>FIGURE 4. Example of time evolution for a geometric Brownian process. We plot the full propagator K given in Eq. (111) for different values of the time $t = 0.2k$, with $k = 1, ..., 10$ (start time in red and end time in blue). We adopted the parameters, $G_0 = 1$, $y = 5$, $t_0 = 0$, and $\alpha = 1/2$ (in arbitrary units).</figcaption></figure>

which turns out to be

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{1}{4(t-t_0)}\left(\int_y^x \frac{d\eta}{g(\eta)}\right)^2\right]}{g(x)\sqrt{4\pi(t-t_0)}}. \quad (110)$$

The Wiener process is of course retrieved when $g(x)$ is a constant.

The fact that in this specific case we were able to find a rather general solution should not be too surprising. In fact, the stochastic differential equation is interpreted by Fisk-Stratonovich integration and thus all the rules of mathematical analysis remain unchanged and applicable as if the equation were not stochastic. Based on this observation, we see that the simple transformation law $z(x) = \int_a^x d\eta/g(\eta)$ (for some real constant $a$) is able to convert the starting equation into $dz/dt = \xi(t)$, which describes a simple Wiener process. This alternative approach leads indeed to the same result stated in Eq. (110). An interesting investigation on this type of equation can be found in Ref. [66].

This problem is useful to further explore the behavior of several systems with heterogeneous diffusion. For example, if we consider the geometric Brownian motion case from Section III with $g(x) = G_0|x|$ but without drift, we simply obtain

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{\ln^2\left|\frac{x}{y}\right|}{4G_0^2(t-t_0)}\right]}{G_0|x|\sqrt{4\pi(t-t_0)}} \mathbf{1}(x), \quad y > 0, \quad (111)$$

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{\ln^2\left|\frac{x}{y}\right|}{4G_0^2(t-t_0)}\right]}{G_0|x|\sqrt{4\pi(t-t_0)}} \mathbf{1}(-x), \quad y < 0, \quad (112)$$

in full agreement with Eqs.(88) and (89), with $H_0 = 0$ and

$\alpha = 1/2$. The geometric Brownian motion without drift, in the case of Fisk-Stratonovich interpretation, can thus be studied directly with this simple approach, without developing the full procedure as in the previous Section III B. An example can be found in Fig. 4, where we plot the time evolution of the propagator given in Eq. (111), with a positive initial condition.

A further interesting case concerns the diffusive behavior described by the power law $g(x) = G_0|x|^{\beta}$, where we assume that $\beta < 1$ for reasons that will be clarified shortly. We have to study the integral $\int_y^x d\eta/|\eta|^{\beta}$. It is rather simple to verify the validity of the indefinite integral $\int d\eta/|\eta|^{\beta} = \operatorname{sgn}(\eta)|\eta|^{1-\beta}/(1-\beta) + C$ for some real constant $C$. Hence, the hypothesis $\beta < 1$ allows us to state that the improper integral is convergent in the neighborhood of $\eta = 0$ (indeed, $1-\beta > 0$), and we can write

$$\begin{aligned} \int_y^x \frac{d\eta}{|\eta|^{\beta}} &= \operatorname{sgn}(x) \frac{|x|^{1-\beta}}{1-\beta} - \operatorname{sgn}(y) \frac{|y|^{1-\beta}}{1-\beta} \\ &= \frac{|x|^{2-\beta}}{x(1-\beta)} - \frac{|y|^{2-\beta}}{y(1-\beta)}, \end{aligned} \quad (113)$$

which is valid for any real numbers $x$ and $y$. Unlike the case of geometric Brownian motion, we find here a propagator defined on the entire real axis (with $\beta < 1$)

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{1}{4G_0^2(t-t_0)(1-\beta)^2}\left(\frac{|x|^{2-\beta}}{x} - \frac{|y|^{2-\beta}}{y}\right)^2\right]}{G_0|x|^{\beta}\sqrt{4\pi(t-t_0)}}. \quad (114)$$

This result represents a generalization of the expression obtained in Ref. [5] for $y = 0$, which reads as

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{1}{4G_0^2(t-t_0)(1-\beta)^2}|x|^{2(1-\beta)}\right]}{G_0|x|^{\beta}\sqrt{4\pi(t-t_0)}}. \quad (115)$$

This process is called subdiffusive if $\beta < 0$ and superdiffusive if $0 < \beta < 1$. This terminology derives from the calculation of mean-squared displacement from Eq. (115), resulting in $\mathbb{E}(x^2(t)) \sim t^{1/(1-\beta)}$, where we identify the exponent $e = 1/(1-\beta)$. Depending on the value of the anomalous diffusion exponent $e$ we can thus distinguish subdiffusion ($0 < e < 1$), superdiffusion ($e > 1$), standard Brownian motion ($e = 1$), ballistic wave-like motion ($e = 2$) and geometric Brownian motion ($e \to \infty$). This analysis shows that the anomalous diffusion observed in several systems can be explained by means of a heterogeneous diffusion, for example of the power-law type. Examples of asymmetric ($y \neq 0$) and symmetric ($y = 0$) subdiffusion and superdiffusion processes can be found in Figs.5 and 6, respectively. In these cases the mean-squared displacement exponent assumes the values $e = 2/3$ (subdiffusive behavior) and $e = 2$ (superdiffusive behavior).

This procedure also works, of course, for the case of heterogeneous exponential diffusion, already discussed in Section III D. The application of Eq. (110) leads directly to the propa-
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Example of time evolution for a geometric Brownian process. We plot the full propagator K given in Eq. (111) for different values of the time $t = 0.2k$, with $k = 1, ..., 10$ (start time in red and end time in blue). We adopted the parameters, $G_0 = 1$, $y = 5$, $t_0 = 0$, and $\alpha = 1/2$ (in arbitrary units).</figcaption></figure>

which turns out to be

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{1}{4(t-t_0)}\left(\int_y^x \frac{d\eta}{g(\eta)}\right)^2\right]}{g(x)\sqrt{4\pi(t-t_0)}}. \quad (110)$$

The Wiener process is of course retrieved when $g(x)$ is a constant.

The fact that in this specific case we were able to find a rather general solution should not be too surprising. In fact, the stochastic differential equation is interpreted by Fisk-Stratonovich integration and thus all the rules of mathematical analysis remain unchanged and applicable as if the equation were not stochastic. Based on this observation, we see that the simple transformation law $z(x) = \int_a^x d\eta/g(\eta)$ (for some real constant $a$) is able to convert the starting equation into $dz/dt = \xi(t)$, which describes a simple Wiener process. This alternative approach leads indeed to the same result stated in Eq. (110). An interesting investigation on this type of equation can be found in Ref. [66].

This problem is useful to further explore the behavior of several systems with heterogeneous diffusion. For example, if we consider the geometric Brownian motion case from Section III with $g(x) = G_0|x|$ but without drift, we simply obtain

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{\ln^2\left|\frac{x}{y}\right|}{4G_0^2(t-t_0)}\right]}{G_0|x|\sqrt{4\pi(t-t_0)}} \mathbf{1}(x), \quad y > 0, \quad (111)$$

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{\ln^2\left|\frac{x}{y}\right|}{4G_0^2(t-t_0)}\right]}{G_0|x|\sqrt{4\pi(t-t_0)}} \mathbf{1}(-x), \quad y < 0, (112)$$

in full agreement with Eqs.(88) and (89), with $H_0 = 0$ and

$\alpha = 1/2$. The geometric Brownian motion without drift, in the case of Fisk-Stratonovich interpretation, can thus be studied directly with this simple approach, without developing the full procedure as in the previous Section III B. An example can be found in Fig. 4, where we plot the time evolution of the propagator given in Eq. (111), with a positive initial condition.

A further interesting case concerns the diffusive behavior described by the power law $g(x) = G_0|x|^{\beta}$, where we assume that $\beta < 1$ for reasons that will be clarified shortly. We have to study the integral $\int_y^x d\eta/|\eta|^{\beta}$. It is rather simple to verify the validity of the indefinite integral $\int d\eta/|\eta|^{\beta} = \operatorname{sgn}(\eta)|\eta|^{1-\beta}/(1-\beta) + C$ for some real constant $C$. Hence, the hypothesis $\beta < 1$ allows us to state that the improper integral is convergent in the neighborhood of $\eta = 0$ (indeed, $1-\beta > 0$), and we can write

$$\begin{aligned} \int_y^x \frac{d\eta}{|\eta|^{\beta}} &= \operatorname{sgn}(x) \frac{|x|^{1-\beta}}{1-\beta} - \operatorname{sgn}(y) \frac{|y|^{1-\beta}}{1-\beta} \\ &= \frac{|x|^{2-\beta}}{x(1-\beta)} - \frac{|y|^{2-\beta}}{y(1-\beta)}, \end{aligned} \quad (113)$$

which is valid for any real numbers $x$ and $y$. Unlike the case of geometric Brownian motion, we find here a propagator defined on the entire real axis (with $\beta < 1$)

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{1}{4G_0^2(t-t_0)(1-\beta)^2}\left(\frac{|x|^{2-\beta}}{x} - \frac{|y|^{2-\beta}}{y}\right)^2\right]}{G_0|x|^{\beta}\sqrt{4\pi(t-t_0)}}. \quad (114)$$

This result represents a generalization of the expression obtained in Ref. [5] for $y = 0$, which reads as

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{1}{4G_0^2(t-t_0)(1-\beta)^2}\left|x^{2(1-\beta)}\right|\right]}{G_0|x|^{\beta}\sqrt{4\pi(t-t_0)}}. \quad (115)$$

This process is called subdiffusive if $\beta < 0$ and superdiffusive if $0 < \beta < 1$. This terminology derives from the calculation of mean-squared displacement from Eq. (115), resulting in $\mathbb{E}(x^2(t)) \sim t^{1/(1-\beta)}$, where we identify the exponent $e = 1/(1-\beta)$. Depending on the value of the anomalous diffusion exponent $e$ we can thus distinguish subdiffusion ($0 < e < 1$), superdiffusion ($e > 1$), standard Brownian motion ($e = 1$), ballistic wave-like motion ($e = 2$) and geometric Brownian motion ($e \to \infty$). This analysis shows that the anomalous diffusion observed in several systems can be explained by means of a heterogeneous diffusion, for example of the power-law type. Examples of asymmetric ($y \neq 0$) and symmetric ($y = 0$) subdiffusion and superdiffusion processes can be found in Figs.5 and 6, respectively. In these cases the mean-squared displacement exponent assumes the values $e = 2/3$ (subdiffusive behavior) and $e = 2$ (superdiffusive behavior).

This procedure also works, of course, for the case of heterogeneous exponential diffusion, already discussed in Section III D. The application of Eq. (110) leads directly to the propa-
<figure><img src="image_5.png" /><figcaption>FIGURE 5. Asymmetric and symmetric subdiffusion processes. We plot the full propagators K given in Eqs.(114) and (115) for different values of the time $t = 0.5k$, with $k = 1, ..., 10$ (start time in red and end time in blue), and for $\beta = -1/2$ ($e = 2/3$). We adopted the parameters, $G_0 = 1$, $y = 3$ (for the asymmetric subdiffusion, left panel), $y = 0$ (for the symmetric subdiffusion, right panel), $t_0 = 0$, and $\alpha = 1/2$ (in arbitrary units).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 6. Asymmetric and symmetric superdiffusion processes. We plot the full propagators K given in Eqs.(114) and (115) for different values of the time $t = 0.5k$ (start time in red and end time in blue), with $k = 1, ..., 10$, and for $\beta = +1/2$ ($e = 2$). We adopted the parameters, $G_0 = 1$, $y = 3$ (for the asymmetric superdiffusion, left panel), $y = 0$ (for the symmetric superdiffusion, right panel), $t_0 = 0$, and $\alpha = 1/2$ (in arbitrary units).</figcaption></figure>

gator in the form

$$K(x, t; y, t_0) = \frac{\exp\left[-\frac{1}{4G_0^2\gamma^2(t-t_0)}(\exp(-\gamma y) - \exp(-\gamma x))^2\right]}{G_0 \exp(\gamma x)\sqrt{4\pi(t-t_0)}}, \quad (116)$$

in agreement with the result in Ref. [5]. In this case, the mean-squared displacement shows a logarithmic time dependence, which is a characteristic feature of ultraslow processes [5]. The existence of the solution for the case with the exponential heterogeneous diffusion under the Fisk-Stratonovich interpretation is consistent with the result of Section III D, affirming that the development in a power series is possible only if $\alpha = 1/2$. An example can be seen in Fig. 7, where we show

the time behavior of the propagator for the exponential process with $G_0 = 1$ and $\gamma = 1$. We observe that the probability density evolves slowly in the direction where the diffusion coefficient is lower, being driven by the gradient of fluctuations.

Although the analysis developed here concerns the Fisk-Stratonovich interpretation, some generalizations can be mentioned. Indeed, it is interesting to remark that Eq.(104), interpreted with $\alpha = 1/2$, is equivalent to the equation

$$\frac{dx}{dt} = (1 - 2\alpha)g(x)\frac{dg(x)}{dx} + g(x)\xi(t), \quad (117)$$

interpreted with an arbitrary value of $\alpha$. Then the solution given in Eq. (110) is valid also in these cases, that are outside
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Asymmetric and symmetric subdiffusion processes. We plot the full propagators K given in Eqs.(114) and (115) for different values of the time $t = 0.5k$, with $k = 1, ..., 10$ (start time in red and end time in blue), and for $\beta = -1/2$ ($e = 2/3$). We adopted the parameters, $G_0 = 1$, $y = 3$ (for the asymmetric subdiffusion, left panel), $y = 0$ (for the symmetric subdiffusion, right panel), $t_0 = 0$, and $\alpha = 1/2$ (in arbitrary units).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Asymmetric and symmetric superdiffusion processes. We plot the full propagators K given in Eqs.(114) and (115) for different values of the time $t = 0.5k$ (start time in red and end time in blue), with $k = 1, ..., 10$, and for $\beta = +1/2$ ($e = 2$). We adopted the parameters, $G_0 = 1$, $y = 3$ (for the asymmetric superdiffusion, left panel), $y = 0$ (for the symmetric superdiffusion, right panel), $t_0 = 0$, and $\alpha = 1/2$ (in arbitrary units).</figcaption></figure>

gator in the form

$$K(x, t; y, t_0) = \frac{\exp\left[-\frac{1}{4G_0^2\gamma^2(t-t_0)}(\exp(-\gamma y) - \exp(-\gamma x))^2\right]}{G_0 \exp(\gamma x)\sqrt{4\pi(t-t_0)}}, \quad (116)$$

in agreement with the result in Ref. [5]. In this case, the mean-squared displacement shows a logarithmic time dependence, which is a characteristic feature of ultraslow processes [5]. The existence of the solution for the case with the exponential heterogeneous diffusion under the Fisk-Stratonovich interpretation is consistent with the result of Section III D, affirming that the development in a power series is possible only if $\alpha = 1/2$. An example can be seen in Fig. 7, where we show

the time behavior of the propagator for the exponential process with $G_0 = 1$ and $\gamma = 1$. We observe that the probability density evolves slowly in the direction where the diffusion coefficient is lower, being driven by the gradient of fluctuations.

Although the analysis developed here concerns the Fisk-Stratonovich interpretation, some generalizations can be mentioned. Indeed, it is interesting to remark that Eq.(104), interpreted with $\alpha = 1/2$, is equivalent to the equation

$$\frac{dx}{dt} = (1 - 2\alpha)g(x)\frac{dg(x)}{dx} + g(x)\xi(t), \quad (117)$$

interpreted with an arbitrary value of $\alpha$. Then the solution given in Eq. (110) is valid also in these cases, that are outside
<figure><img src="image_7.png" /><figcaption>FIGURE 7. Example of time evolution for an exponential process. We plot the full propagator K given in Eq. (116) for different values of the time $t = 1, ..., 10$ (start time in red and end time in blue). We adopted the parameter values $G_0 = 1$, $\gamma = 1$, $y = 0$, $t_0 = 0$, and $\alpha = 1/2$ (in arbitrary units).</figcaption></figure>

the Fisk-Stratonovich interpretation.

## B. Stochastic processes with drift under the Fisk-Stratonovich interpretation

We now return to stochastic processes with drift, studied under the Fisk-Stratonovich interpretation, and search for functions $h$ and $g$ that are able to significantly simplify the recursive procedure obtained with our approach. In particular, we look for the form of these functions that generate constant coefficients $\Lambda_n$ and $D_n$ (i.e., independent of $x$). From Eq.(63), we see this can be true only if

$$ -f - \ell\Phi + g^2 \frac{\partial\Phi}{\partial x} + g^2\Phi^2 = C, \quad (118) $$

where $C$ is a constant. By using Eqs.(5), (6), and (58) for $f$, $\ell$, and $\Phi$, this condition can be rewritten as

$$ \begin{aligned} & h' + 2(\alpha - 1)h \frac{g'}{g} + \frac{h^2}{2g^2} \\ & = -2C - (2\alpha - 1) \left[ \frac{2\alpha - 1}{2} (g')^2 + gg'' \right], \end{aligned} \quad (119) $$

where $\alpha$ is the stochastic discretization parameter. This relation must be fulfilled by $h$ and $g$ in order to have constant coefficients $\Lambda_n$ and $D_n$. In the particular case with $h = 0$ and $\alpha = 1/2$, we get $C = 0$, and therefore we retrieve the results of the previous Section (all coefficients are zero for $n \ge 1$). We now want to retain the Fisk-Stratonovich interpretation but

also assume that $h \neq 0$. In this situation, Eq.(119) reduces to

$$ h' - h \frac{g'}{g} + \frac{h^2}{2g^2} = -2C, \quad (120) $$

A simple assumption that allows us to verify this condition is given by $h(x) = -mg(x)$, where $m$ is a real parameter. It means that we are studying the stochastic differential equation

$$ \frac{dx}{dt} = -mg(x) + g(x)\xi(t), \quad (121) $$

where the drift term is proportional to the diffusion term. In this particular case we get

$$ f(x) = -mg' - (g')^2 - gg'', \quad (122) $$

$$ \ell(x) = -mg - 3gg', \quad (123) $$

$$ \Phi(x) = -\frac{1}{2} \frac{m}{g} - \frac{g'}{g}. \quad (124) $$

The short-time propagator is given by the expression

$$ K_0(x,t;y,t_0) = \frac{\exp\left[-\frac{1}{4(t-t_0)}\left(\int_y^x \frac{d\eta}{g(\eta)}\right)^2\right]}{g(y)\sqrt{4\pi(t-t_0)}}, \quad (125) $$

The first correction term $F_0$ can be calculated as follows

$$ F_0 = \exp\left[\int_y^x \Phi(\eta) d\eta\right] = \frac{g(y)}{g(x)} \exp\left[-\frac{m}{2} \int_y^x \frac{d\eta}{g(\eta)}\right], \quad (126) $$

and the other coefficients must be determined through the iteration stated in Eqs.(63) and (65). In this case, Eq.(118) is fulfilled with a constant $C = -m^2/4$, and therefore we have from Eq.(63)

$$ \Lambda_n = -\frac{1}{4}m^2 D_{n-1}. \quad (127) $$

When this relation is substituted into Eq.(65), the following simple recursion is found

$$ D_n = -\frac{1}{4n}m^2 D_{n-1}, \quad (128) $$

which leads to the result

$$ D_n = (-1)^n \frac{m^{2n}}{4^n n!} \quad (129) $$

that complies with the condition $D_0 = 1$. The correction function $F$ can be summed as follows

$$ F = F_0 \sum_{n=0}^{+\infty} \frac{1}{n!} \left[ -\frac{m^2 (t-t_0)}{4} \right]^n = F_0 \exp \left[ -\frac{m^2}{4} (t-t_0) \right], \quad (130) $$
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Example of time evolution for an exponential process. We plot the full propagator K given in Eq. (116) for different values of the time $t = 1, ..., 10$ (start time in red and end time in blue). We adopted the parameter values $G_0 = 1$, $\gamma = 1$, $y = 0$, $t_0 = 0$, and $\alpha = 1/2$ (in arbitrary units).</figcaption></figure>

the Fisk-Stratonovich interpretation.

## B. Stochastic processes with drift under the Fisk-Stratonovich interpretation

We now return to stochastic processes with drift, studied under the Fisk-Stratonovich interpretation, and search for functions $h$ and $g$ that are able to significantly simplify the recursive procedure obtained with our approach. In particular, we look for the form of these functions that generate constant coefficients $\Lambda_n$ and $D_n$ (i.e., independent of $x$). From Eq.(63), we see this can be true only if

$$ -f - \ell\Phi + g^2 \frac{\partial\Phi}{\partial x} + g^2\Phi^2 = C, \quad (118) $$

where $C$ is a constant. By using Eqs.(5), (6), and (58) for $f$, $\ell$, and $\Phi$, this condition can be rewritten as

$$ \begin{aligned} & h' + 2(\alpha - 1)h \frac{g'}{g} + \frac{h^2}{2g^2} \\ & = -2C - (2\alpha - 1) \left[ \frac{2\alpha - 1}{2} (g')^2 + gg'' \right], \end{aligned} \quad (119) $$

where $\alpha$ is the stochastic discretization parameter. This relation must be fulfilled by $h$ and $g$ in order to have constant coefficients $\Lambda_n$ and $D_n$. In the particular case with $h = 0$ and $\alpha = 1/2$, we get $C = 0$, and therefore we retrieve the results of the previous Section (all coefficients are zero for $n \ge 1$). We now want to retain the Fisk-Stratonovich interpretation but

also assume that $h \neq 0$. In this situation, Eq.(119) reduces to

$$ h' - h \frac{g'}{g} + \frac{h^2}{2g^2} = -2C, \quad (120) $$

A simple assumption that allows us to verify this condition is given by $h(x) = -mg(x)$, where $m$ is a real parameter. It means that we are studying the stochastic differential equation

$$ \frac{dx}{dt} = -mg(x) + g(x)\xi(t), \quad (121) $$

where the drift term is proportional to the diffusion term. In this particular case we get

$$ f(x) = -mg' - (g')^2 - gg'', \quad (122) $$

$$ \ell(x) = -mg - 3gg', \quad (123) $$

$$ \Phi(x) = -\frac{1}{2} \frac{m}{g} - \frac{g'}{g}. \quad (124) $$

The short-time propagator is given by the expression

$$ K_0(x, t; y, t_0) = \frac{\exp\left[-\frac{1}{4(t-t_0)}\left(\int_y^x \frac{d\eta}{g(\eta)}\right)^2\right]}{g(y)\sqrt{4\pi(t-t_0)}}, \quad (125) $$

The first correction term $F_0$ can be calculated as follows

$$ F_0 = \exp\left[\int_y^x \Phi(\eta) d\eta\right] = \frac{g(y)}{g(x)} \exp\left[-\frac{m}{2} \int_y^x \frac{d\eta}{g(\eta)}\right], \quad (126) $$

and the other coefficients must be determined through the iteration stated in Eqs.(63) and (65). In this case, Eq.(118) is fulfilled with a constant $C = -m^2/4$, and therefore we have from Eq.(63)

$$ \Lambda_n = -\frac{1}{4}m^2 D_{n-1}. \quad (127) $$

When this relation is substituted into Eq.(65), the following simple recursion is found

$$ D_n = -\frac{1}{4n}m^2 D_{n-1}, \quad (128) $$

which leads to the result

$$ D_n = (-1)^n \frac{m^{2n}}{4^n n!} \quad (129) $$

that complies with the condition $D_0 = 1$. The correction function $F$ can be summed as follows

$$ F = F_0 \sum_{n=0}^{+\infty} \frac{1}{n!} \left[ -\frac{m^2 (t-t_0)}{4} \right]^n = F_0 \exp \left[ -\frac{m^2}{4} (t-t_0) \right], \quad (130) $$
and the propagator $K = K_0 F$ is finally found in the form

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{1}{4(t-t_0)}\left(\int_y^x \frac{d\eta}{g(\eta)} + m(t-t_0)\right)^2\right]}{g(x)\sqrt{4\pi(t-t_0)}}. \quad (131)$$

This result, as expected, reduce to Eq.(110) when $m = 0$, i.e., when we set the drift to zero. The problem described by Eq.(121), just solved for $\alpha = 1/2$, is useful for studying various stochastic population growth models. For instance, as the diffusion function $g$ varies, we can recover the Verhulst logistic model, the Gompertz model, the Shoener model, the Richards model, and the Smith model [67–71].

Another interesting example can be found by considering $g(x)$ fixed and by searching for the solution of Eq.(120) with respect to $h(x)$. This equation is of Riccati type, which means that it is composed of a linear combination of the derivative of $h$, a linear term in $h$, a term independent of $h$ and finally a quadratic term in $h$. We learned from the previous example that $h(x) = -mg(x)$ is a particular solution of this equation (with $C = -m^2/4$). We can therefore search for a general solution of Eq.(120) in the form $h(x) = -mg(x) + \tau(x)$ (again with $C = -m^2/4$). A simple calculation proves that $\tau(x)$ satisfies the following equation

$$\tau' = \left(\frac{g'}{g} + \frac{m}{g}\right) \tau - \frac{\tau^2}{2g^2}, \quad (132)$$

which is of the Bernoulli type (with a structure similar to the Riccati equation without the term independent of the unknown $\tau$). It means that the knowledge of a particular solution of Eq.(120), in the form $h(x) = -mg(x)$, allowed us to move from a Riccati equation to a Bernoulli equation, thus slightly simpler. A further change of variable, namely $\theta = 1/\tau$, finally reduces the equation to the linear case. Straightforward calculations lead indeed to the equation

$$\theta' = - \left( \frac{g'}{g} + \frac{m}{g} \right) \theta + \frac{1}{2g^2}, \quad (133)$$

which can be solved without additional effort. In fact, by defining a function $p(x)$ such that $p'(x) = 1/g(x)$, i.e. $p$ is any indefinite integral of $1/g$, we have the solution

$$\theta(x) = \frac{S}{g(x)} e^{-mp(x)} + \frac{1}{2g(x)} e^{-mp(x)} \int^{x} \frac{1}{g(\eta)} e^{mp(\eta)} d\eta, \quad (134)$$

where $S$ is an arbitrary integration constant. We can conclude that the general solution of Eq.(120) is given by $h(x) = -mg(x) + 1/\theta(x)$. As a result of the construction developed, the Fokker-Planck equation generated by $h$ and $g$ definitely has constant coefficients $D_n$ and thus its propagator can be found without difficulty. This defines a new class of Langevin equations for which the exact solution can be calculated.

Let's study a specific example in detail concerning a constant diffusion term $g = G_0$. This choice leads to a drift term $h(x) = -mG_0 + 1/\theta(x) = -H_0 + 1/\theta(x)$, where $\theta(x)$

can be obtained through Eqs.(133) and (134). More explicitly

$$h(x) = H_0 \frac{1 - 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)}{1 + 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)}, \quad (135)$$

where $S$ is arbitrary. We are therefore studying the nonlinear Langevin equation with a sigmoidal drift

$$\frac{dx}{dt} = H_0 \frac{1 - 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)}{1 + 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)} + G_0 \xi(t), \quad (136)$$

interpreted through $\alpha = 1/2$. For this equation, the function $\Phi(x)$ can be written as

$$\Phi(x) = \frac{h(x)}{2G_0^2} = \frac{H_0}{2G_0^2} \frac{1 - 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)}{1 + 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)}, \quad (137)$$

and therefore

$$F_0 = \exp\left[\int_y^x \Phi(\eta) d\eta\right] = \frac{\exp\left(\frac{H_0}{2G_0^2}x\right) + 2H_0 S \exp\left(-\frac{H_0}{2G_0^2}x\right)}{\exp\left(\frac{H_0}{2G_0^2}y\right) + 2H_0 S \exp\left(-\frac{H_0}{2G_0^2}y\right)}. \quad (138)$$

Moreover, the constant $C$ is given by $C = -m^2/4 = -H_0^2/(4G_0^2)$. Hence, we have

$$D_n = (-1)^n \frac{m^{2n}}{4^n n!} = (-1)^n \frac{1}{4^n n!} \left(\frac{H_0}{G_0}\right)^{2n}. \quad (139)$$

We can calculate the correction function $F$

$$F = F_0 \exp \left[ -\frac{1}{4} \left( \frac{H_0}{G_0} \right)^2 (t-t_0) \right], \quad (140)$$

and we finally have the propagator $K = K_0 F$ in explicit form

$$K(x,t;y,t_0) = \frac{1}{\sqrt{4\pi G_0^2 (t-t_0)}} \exp\left[-\frac{(x-y)^2}{4G_0^2 (t-t_0)}\right] \\ \times \frac{\exp\left(\frac{H_0}{2G_0^2}x\right) + 2H_0 S \exp\left(-\frac{H_0}{2G_0^2}x\right)}{\exp\left(\frac{H_0}{2G_0^2}y\right) + 2H_0 S \exp\left(-\frac{H_0}{2G_0^2}y\right)} \\ \times \exp\left[-\frac{1}{4}\left(\frac{H_0}{G_0}\right)^2 (t-t_0)\right]. \quad (141)$$

It can be verified that this expression is indeed the solution of the Fokker-Planck equation associated to Eq.(136). A numerical example can be found in Fig. 8, where we show the sigmoidal drift with the propagator evolution for two different initial conditions, corresponding to the positive and negative regions of the sigmoidal function. We observe a regressive diffusion in the first case and a progressive diffusion in the second one, as
and the propagator $K = K_0 F$ is finally found in the form

$$K(x,t;y,t_0) = \frac{\exp\left[-\frac{1}{4(t-t_0)}\left(\int_y^x \frac{d\eta}{g(\eta)} + m(t-t_0)\right)^2\right]}{g(x)\sqrt{4\pi(t-t_0)}}. \quad (131)$$

This result, as expected, reduce to Eq.(110) when $m = 0$, i.e., when we set the drift to zero. The problem described by Eq.(121), just solved for $\alpha = 1/2$, is useful for studying various stochastic population growth models. For instance, as the diffusion function $g$ varies, we can recover the Verhulst logistic model, the Gompertz model, the Shoener model, the Richards model, and the Smith model [67–71].

Another interesting example can be found by considering $g(x)$ fixed and by searching for the solution of Eq.(120) with respect to $h(x)$. This equation is of Riccati type, which means that it is composed of a linear combination of the derivative of $h$, a linear term in $h$, a term independent of $h$ and finally a quadratic term in $h$. We learned from the previous example that $h(x) = -mg(x)$ is a particular solution of this equation (with $C = -m^2/4$). We can therefore search for a general solution of Eq.(120) in the form $h(x) = -mg(x) + \tau(x)$ (again with $C = -m^2/4$). A simple calculation proves that $\tau(x)$ satisfies the following equation

$$\tau' = \left(\frac{g'}{g} + \frac{m}{g}\right) \tau - \frac{\tau^2}{2g^2}, \quad (132)$$

which is of the Bernoulli type (with a structure similar to the Riccati equation without the term independent of the unknown $\tau$). It means that the knowledge of a particular solution of Eq.(120), in the form $h(x) = -mg(x)$, allowed us to move from a Riccati equation to a Bernoulli equation, thus slightly simpler. A further change of variable, namely $\theta = 1/\tau$, finally reduces the equation to the linear case. Straightforward calculations lead indeed to the equation

$$\theta' = -\left(\frac{g'}{g} + \frac{m}{g}\right) \theta + \frac{1}{2g^2}, \quad (133)$$

which can be solved without additional effort. In fact, by defining a function $p(x)$ such that $p'(x) = 1/g(x)$, i.e. $p$ is any indefinite integral of $1/g$, we have the solution

$$\theta(x) = \frac{S}{g(x)} e^{-mp(x)} + \frac{1}{2g(x)} e^{-mp(x)} \int^{x} \frac{1}{g(\eta)} e^{mp(\eta)} d\eta, \quad (134)$$

where $S$ is an arbitrary integration constant. We can conclude that the general solution of Eq.(120) is given by $h(x) = -mg(x) + 1/\theta(x)$. As a result of the construction developed, the Fokker-Planck equation generated by $h$ and $g$ definitely has constant coefficients $D_n$ and thus its propagator can be found without difficulty. This defines a new class of Langevin equations for which the exact solution can be calculated.

Let's study a specific example in detail concerning a constant diffusion term $g = G_0$. This choice leads to a drift term $h(x) = -mG_0 + 1/\theta(x) = -H_0 + 1/\theta(x)$, where $\theta(x)$

can be obtained through Eqs.(133) and (134). More explicitly

$$h(x) = H_0 \frac{1 - 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)}{1 + 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)}, \quad (135)$$

where $S$ is arbitrary. We are therefore studying the nonlinear Langevin equation with a sigmoidal drift

$$\frac{dx}{dt} = H_0 \frac{1 - 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)}{1 + 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)} + G_0 \xi(t), \quad (136)$$

interpreted through $\alpha = 1/2$. For this equation, the function $\Phi(x)$ can be written as

$$\Phi(x) = \frac{h(x)}{2G_0^2} = \frac{H_0}{2G_0^2} \frac{1 - 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)}{1 + 2H_0 S \exp\left(-\frac{H_0}{G_0^2}x\right)}, \quad (137)$$

and therefore

$$F_0 = \exp\left[\int_y^x \Phi(\eta) d\eta\right] = \frac{\exp\left(\frac{H_0}{2G_0^2}x\right) + 2H_0 S \exp\left(-\frac{H_0}{2G_0^2}x\right)}{\exp\left(\frac{H_0}{2G_0^2}y\right) + 2H_0 S \exp\left(-\frac{H_0}{2G_0^2}y\right)}. \quad (138)$$

Moreover, the constant $C$ is given by $C = -m^2/4 = -H_0^2/(4G_0^2)$. Hence, we have

$$D_n = (-1)^n \frac{m^{2n}}{4^n n!} = (-1)^n \frac{1}{4^n n!} \left(\frac{H_0}{G_0}\right)^{2n}. \quad (139)$$

We can calculate the correction function $F$

$$F = F_0 \exp\left[-\frac{1}{4}\left(\frac{H_0}{G_0}\right)^2 (t-t_0)\right], \quad (140)$$

and we finally have the propagator $K = K_0 F$ in explicit form

$$K(x,t;y,t_0) = \frac{1}{\sqrt{4\pi G_0^2 (t-t_0)}} \exp\left[-\frac{(x-y)^2}{4G_0^2 (t-t_0)}\right] \\ \times \frac{\exp\left(\frac{H_0}{2G_0^2}x\right) + 2H_0 S \exp\left(-\frac{H_0}{2G_0^2}x\right)}{\exp\left(\frac{H_0}{2G_0^2}y\right) + 2H_0 S \exp\left(-\frac{H_0}{2G_0^2}y\right)} \\ \times \exp\left[-\frac{1}{4}\left(\frac{H_0}{G_0}\right)^2 (t-t_0)\right]. \quad (141)$$

It can be verified that this expression is indeed the solution of the Fokker-Planck equation associated to Eq.(136). A numerical example can be found in Fig. 8, where we show the sigmoidal drift with the propagator evolution for two different initial conditions, corresponding to the positive and negative regions of the sigmoidal function. We observe a regressive diffusion in the first case and a progressive diffusion in the second one, as
<figure><img src="image_8.png" /><figcaption>FIGURE 8. Stochastic process with a sigmoidal drift. We plot the full propagators K given in Eq.(141) for different values of the time $t = 0.1k$, with $k = 1, ..., 20$ (start time in red and end time in blue). We also plot the sigmoidal drift $h(x)$ (black lines). We adopted the parameter values, $G_0 = 1$, $H_0 = 1$, $t_0 = 0$, $S = 10$, and $\alpha = 1/2$ (in arbitrary units). We observe a regressive diffusion in left panel with $y = 0$, and a progressive diffusion in the right panel with $y = 5$.</figcaption></figure>

expected. We also observe that the rate of displacement of the mean value of the variable $x$ converges to a well-determined value for long times because of the horizontal asymptotes of the sigmoidal function. Moreover, the variance is an increasing function due to the noise fluctuations.

Another example deals with the heterogeneous exponential diffusion of the form $g(x) = G_0 \exp(\gamma x)$. We consider again $\alpha = 1/2$, and we search for the function $h$ given by $h(x) = -mg(x) + 1/\theta(x)$, as before. We define $H_0 = mG_0$, and we obtain $p(x) = -\exp(-\gamma x)/(G_0\gamma)$, such that $p'(x) = 1/g(x)$. The application of previous procedure delivers

$$h(x) = H_0 \exp(\gamma x) \frac{1 - 2 \frac{H_0}{G_0} S \exp \left[ -\frac{H_0}{G_0^2} p(x) \right]}{1 + 2 \frac{H_0}{G_0} S \exp \left[ -\frac{H_0}{G_0^2} p(x) \right]} \\ = H_0 \exp(\gamma x) \frac{1 - 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma x) \right]}{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma x) \right]}, \quad (142)$$

where $S$ is an arbitrary constant. It means that we study the nonlinear Langevin equation with heterogeneous diffusion

$$\frac{dx}{dt} = H_0 \exp(\gamma x) \frac{1 - 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma x) \right]}{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma x) \right]} + G_0 \exp(\gamma x) \xi(t), \quad (143)$$

under the Fisk-Stratonovich interpretation. With straightforward calculation we find

$$F_0 = \frac{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2 \gamma} \exp(-\gamma x) \right]}{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2 \gamma} \exp(-\gamma y) \right]} \times \exp \left\{ -\gamma(x-y) - \frac{H_0}{2G_0^2 \gamma} \left[ \exp(-\gamma x) - \exp(-\gamma y) \right] \right\}, \quad (144)$$

and the correction function $F$ appears to be

$$F = F_0 \exp \left[ -\frac{1}{4} \left( \frac{H_0}{G_0} \right)^2 (t - t_0) \right]. \quad (145)$$

To conclude, we can obtain the short time propagator as

$$K_0 = \frac{\exp \left\{ - \frac{[\exp(-\gamma x) - \exp(-\gamma y)]^2}{4G_0^2 \gamma^2 (t-t_0)} \right\}}{\sqrt{4\pi G_0^2 \exp(2\gamma y)(t-t_0)}}, \quad (146)$$

and the complete propagator in the form

$$K = \frac{\exp \left\{ - \frac{[\exp(-\gamma x) - \exp(-\gamma y)]^2}{4G_0^2 \gamma^2 (t-t_0)} \right\}}{\sqrt{4\pi G_0^2 \exp(2\gamma y)(t-t_0)}} \exp \left[ -\frac{1}{4} \left( \frac{H_0}{G_0} \right)^2 (t-t_0) \right] \\ \times \frac{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma x) \right]}{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma y) \right]} \\ \times \exp \left\{ -\gamma(x-y) - \frac{H_0}{2G_0^2 \gamma} \left[ \exp(-\gamma x) - \exp(-\gamma y) \right] \right\}. \quad (147)$$

This is the final result concerning the exponential diffusion with a specific drift and it can be checked by direct substitution calculation.
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 8. Stochastic process with a sigmoidal drift. We plot the full propagators K given in Eq.(141) for different values of the time t = 0.1k, with k = 1, ..., 20 (start time in red and end time in blue). We also plot the sigmoidal drift h(x) (black lines). We adopted the parameter values, G₀ = 1, H₀ = 1, t₀ = 0, S = 10, and α = 1/2 (in arbitrary units). We observe a regressive diffusion in left panel with y = 0, and a progressive diffusion in the right panel with y = 5.</figcaption></figure>

expected. We also observe that the rate of displacement of the mean value of the variable x converges to a well-determined value for long times because of the horizontal asymptotes of the sigmoidal function. Moreover, the variance is an increasing function due to the noise fluctuations.

Another example deals with the heterogeneous exponential diffusion of the form $g(x) = G_0 \exp(\gamma x)$. We consider again $\alpha = 1/2$, and we search for the function $h$ given by $h(x) = -mg(x) + 1/\theta(x)$, as before. We define $H_0 = mG_0$, and we obtain $p(x) = -\exp(-\gamma x)/(G_0\gamma)$, such that $p'(x) = 1/g(x)$. The application of previous procedure delivers

$$h(x) = H_0 \exp(\gamma x) \frac{1 - 2 \frac{H_0}{G_0} S \exp \left[ -\frac{H_0}{G_0^2} p(x) \right]}{1 + 2 \frac{H_0}{G_0} S \exp \left[ -\frac{H_0}{G_0^2} p(x) \right]} \\ = H_0 \exp(\gamma x) \frac{1 - 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma x) \right]}{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma x) \right]}, \quad (142)$$

where S is an arbitrary constant. It means that we study the nonlinear Langevin equation with heterogeneous diffusion

$$\frac{dx}{dt} = H_0 \exp(\gamma x) \frac{1 - 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma x) \right]}{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma x) \right]} + G_0 \exp(\gamma x) \xi(t), \quad (143)$$

under the Fisk-Stratonovich interpretation. With straightforward calculation we find

$$F_0 = \frac{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma x) \right]}{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma y) \right]} \times \exp \left\{ -\gamma(x-y) - \frac{H_0}{2G_0^2 \gamma} \left[ \exp(-\gamma x) - \exp(-\gamma y) \right] \right\}, \quad (144)$$

and the correction function F appears to be

$$F = F_0 \exp \left[ -\frac{1}{4} \left( \frac{H_0}{G_0} \right)^2 (t - t_0) \right]. \quad (145)$$

To conclude, we can obtain the short time propagator as

$$K_0 = \frac{\exp \left\{ - \frac{[\exp(-\gamma x) - \exp(-\gamma y)]^2}{4 G_0^2 \gamma^2 (t - t_0)} \right\}}{\sqrt{4 \pi G_0^2 \exp(2\gamma y) (t - t_0)}}, \quad (146)$$

and the complete propagator in the form

$$K = \frac{\exp \left\{ - \frac{[\exp(-\gamma x) - \exp(-\gamma y)]^2}{4 G_0^2 \gamma^2 (t-t_0)} \right\}}{\sqrt{4 \pi G_0^2 \exp(2\gamma y) (t-t_0)}} \exp \left[ -\frac{1}{4} \left( \frac{H_0}{G_0} \right)^2 (t-t_0) \right] \\ \times \frac{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma x) \right]}{1 + 2 \frac{H_0}{G_0} S \exp \left[ \frac{H_0}{G_0^2} \gamma \exp(-\gamma y) \right]} \\ \times \exp \left\{ -\gamma(x-y) - \frac{H_0}{2G_0^2 \gamma} \left[ \exp(-\gamma x) - \exp(-\gamma y) \right] \right\}. \quad (147)$$

This is the final result concerning the exponential diffusion with a specific drift and it can be checked by direct substitution calculation.
tion into the Fokker-Planck equation.

This procedure, once fixed an arbitrary heterogeneous diffusion described by $g(x)$, enables us to find a function $h(x)$ that allows solving the associated Fokker-Planck equation with some degree of ease, as shown in the examples above. This can be useful for constructing stochastic systems whose solutions are known, with a variety of practical applications.

## V. CONCLUSIONS

In this paper we have derived a short-time expansion for the probability distribution of a Fokker-Planck equation. We have shown that its kernel or propagator can be expressed as the product of a singular and a regular term, the latter of which allows an expansion in a Taylor series, whose coefficients can be obtained from the recursive solutions of a system of ordinary differential equations. Their solutions simplify considerably in the case of time-independent drift and diffusion coefficients. Indeed, with fully time-dependent drift and diffusion terms, each Taylor coefficient of the expansion depends on all the previous coefficients. Differently, with time-independent drift and diffusion terms, each coefficient only depends on the previous one. This reduction in complexity, allows the procedure to be applied to cases of practical interest much more easily. Moreover, it can be applied manually in simple cases of theoretically interest and can be implemented in automatic symbolic computing environments in more complex cases of application interest.

Subsequently we have put our newly developed machinery to work on known examples of stochastic processes, highlighting specific aspects of the implementation of our approach. Specifically, we have shown the application of our method to the simple case of Gaussian processes, to geometrical Brownian motions leading to the so-called log-normal distribution, and to a particular stochastic process used to describe chromatin remodeling phenomena. In particular, we see that the expansion cannot be guaranteed to work independently of the discretization parameter $\alpha$, which we illustrated for the case of an exponential diffusion coefficient. We also showed that in the cases in which the expansion works and leads to a convergent expression for the time-dependent propagator - hence also at long times - only a few of the expansion coefficients are needed in order to already obtain a quantitatively reliable result.

Finally, we have taken advantage of our recursion system for the determination of the expansion coefficients to allow us to find novel drift and diffusion terms for which the procedure yields exact results for the propagator. This feature is of interest to be further exploited to define novel stochastic processes of interest in future applications. Within the Fisk-Stratonovich interpretation, this approach allowed us to solve the stochastic differential equation without drift term. This result provides a unified picture for describing the geometric Brownian motion, the power law diffusion, and the exponential diffusion. These different heterogeneous diffusion schemes are well adapted to describe the anomalous diffusion observed in several physical and biological systems. A further generalization concerns

the Fisk-Stratonovich equation with drift. Here we have proposed a technique that allows us to identify pairs of drift and diffusion terms for which it is easy to determine the exact expression of the propagator. This is useful for exploring new forms of stochastic equations and processes with possible various applications. From this point of view, an interesting future perspective concerns the application of this approach to problems with drift or diffusion terms that are periodic in space. This can be useful to study diffusion processes in stratified or multilayered systems and to analyze the behaviors of reaction-diffusion equations with application to chemical processes and other physical instabilities and transformations.

## VI. APPENDIX

In this appendix we list, for completeness, the coefficients of the development for the model discussed in Section III C, with the results shown in Fig. 3.

A first set of parameters is given by $k = 2$, $\beta = 1$, $\delta = 1/10$, and $A = 1$ (in arbitrary units) and the results have been obtained as

$$D_0 = 1, \qquad (148)$$

$$D_1 = 2 + \frac{1}{10}x - \frac{67}{50}x^2 - \frac{1099}{3000}x^3 - \frac{409}{10000}x^4 - \frac{91}{45000}x^5 + O(x^6), \qquad (149)$$

$$D_2 = \frac{13}{20} + \frac{1}{150}x - \frac{19969}{7500}x^2 - \frac{3251}{3750}x^3 + \frac{1363801}{1750000}x^4 + \frac{50689969}{105000000}x^5 + O(x^6), \qquad (150)$$

$$D_3 = -\frac{5033}{3750} - \frac{11641}{25000}x + \frac{1521223}{2625000}x^2 + \frac{4000063}{13125000}x^3 + \frac{1149745319}{630000000}x^4 + \frac{14432096843}{13500000000}x^5 + O(x^6), \qquad (151)$$

$$D_4 = -\frac{3759293}{10500000} - \frac{16963}{750000}x + \frac{1454450549}{315000000}x^2 + \frac{11140823723}{4725000000}x^3 - \frac{57641437687}{630000000000}x^4 - \frac{458475394523}{506250000000}x^5 + O(x^6). \qquad (152)$$

A second set of parameters is given by $k = 2$, $\beta = 1$, $\delta = 1/10$, and $A = 2$ (in arbitrary units) and we obtained

$$D_0 = 1, \qquad (153)$$
tion into the Fokker-Planck equation.

This procedure, once fixed an arbitrary heterogeneous diffusion described by $g(x)$, enables us to find a function $h(x)$ that allows solving the associated Fokker-Planck equation with some degree of ease, as shown in the examples above. This can be useful for constructing stochastic systems whose solutions are known, with a variety of practical applications.

## V. CONCLUSIONS

In this paper we have derived a short-time expansion for the probability distribution of a Fokker-Planck equation. We have shown that its kernel or propagator can be expressed as the product of a singular and a regular term, the latter of which allows an expansion in a Taylor series, whose coefficients can be obtained from the recursive solutions of a system of ordinary differential equations. Their solutions simplify considerably in the case of time-independent drift and diffusion coefficients. Indeed, with fully time-dependent drift and diffusion terms, each Taylor coefficient of the expansion depends on all the previous coefficients. Differently, with time-independent drift and diffusion terms, each coefficient only depends on the previous one. This reduction in complexity, allows the procedure to be applied to cases of practical interest much more easily. Moreover, it can be applied manually in simple cases of theoretically interest and can be implemented in automatic symbolic computing environments in more complex cases of application interest.

Subsequently we have put our newly developed machinery to work on known examples of stochastic processes, highlighting specific aspects of the implementation of our approach. Specifically, we have shown the application of our method to the simple case of Gaussian processes, to geometrical Brownian motions leading to the so-called log-normal distribution, and to a particular stochastic process used to describe chromatin remodeling phenomena. In particular, we see that the expansion cannot be guaranteed to work independently of the discretization parameter $\alpha$, which we illustrated for the case of an exponential diffusion coefficient. We also showed that in the cases in which the expansion works and leads to a convergent expression for the time-dependent propagator - hence also at long times - only a few of the expansion coefficients are needed in order to already obtain a quantitatively reliable result.

Finally, we have taken advantage of our recursion system for the determination of the expansion coefficients to allow us to find novel drift and diffusion terms for which the procedure yields exact results for the propagator. This feature is of interest to be further exploited to define novel stochastic processes of interest in future applications. Within the Fisk-Stratonovich interpretation, this approach allowed us to solve the stochastic differential equation without drift term. This result provides a unified picture for describing the geometric Brownian motion, the power law diffusion, and the exponential diffusion. These different heterogeneous diffusion schemes are well adapted to describe the anomalous diffusion observed in several physical and biological systems. A further generalization concerns

the Fisk-Stratonovich equation with drift. Here we have proposed a technique that allows us to identify pairs of drift and diffusion terms for which it is easy to determine the exact expression of the propagator. This is useful for exploring new forms of stochastic equations and processes with possible various applications. From this point of view, an interesting future perspective concerns the application of this approach to problems with drift or diffusion terms that are periodic in space. This can be useful to study diffusion processes in stratified or multilayered systems and to analyze the behaviors of reaction-diffusion equations with application to chemical processes and other physical instabilities and transformations.

## VI. APPENDIX

In this appendix we list, for completeness, the coefficients of the development for the model discussed in Section III C, with the results shown in Fig. 3.

A first set of parameters is given by $k = 2$, $\beta = 1$, $\delta = 1/10$, and $A = 1$ (in arbitrary units) and the results have been obtained as

$$D_0 = 1, \quad (148)$$

$$D_1 = 2 + \frac{1}{10}x - \frac{67}{50}x^2 - \frac{1099}{3000}x^3 - \frac{409}{10000}x^4 - \frac{91}{45000}x^5 + O(x^6), \quad (149)$$

$$D_2 = \frac{13}{20} + \frac{1}{150}x - \frac{19969}{7500}x^2 - \frac{3251}{3750}x^3 + \frac{1363801}{1750000}x^4 + \frac{50689969}{105000000}x^5 + O(x^6), \quad (150)$$

$$D_3 = -\frac{5033}{3750} - \frac{11641}{25000}x + \frac{1521223}{2625000}x^2 + \frac{4000063}{13125000}x^3 + \frac{1149745319}{630000000}x^4 + \frac{14432096843}{13500000000}x^5 + O(x^6), \quad (151)$$

$$D_4 = -\frac{3759293}{10500000} - \frac{16963}{750000}x + \frac{1454450549}{315000000}x^2 + \frac{11140823723}{4725000000}x^3 - \frac{57641437687}{630000000000}x^4 - \frac{458475394523}{506250000000}x^5 + O(x^6). \quad (152)$$

A second set of parameters is given by $k = 2$, $\beta = 1$, $\delta = 1/10$, and $A = 2$ (in arbitrary units) and we obtained

$$D_0 = 1, \quad (153)$$
$$D_1 = 2 + \frac{1}{10}x - \frac{17}{50}x^2 - \frac{137}{1500}x^3 - \frac{307}{30000}x^4 - \frac{91}{180000}x^5 + O(x^6), \quad (154)$$

$$D_2 = \frac{3}{5} + \frac{2}{75}x - \frac{10007}{15000}x^2 - \frac{1621}{7500}x^3 + \frac{24527}{875000}x^4 + \frac{1524793}{52500000}x^5 + O(x^6), \quad (155)$$

$$D_3 = -\frac{2599}{1875} - \frac{1407}{3125}x + \frac{57521}{328125}x^2 + \frac{1992491}{26250000}x^3 + \frac{9657743}{78750000}x^4 + \frac{795208357}{11812500000}x^5 + O(x^6), \quad (156)$$

$$D_4 = -\frac{165647}{656250} - \frac{7079}{93750}x + \frac{45878477}{39375000}x^2 + \frac{347956369}{590625000}x^3 + \frac{741347219}{47250000000}x^4 - \frac{45979054151}{8859375000000}x^5 + O(x^6), \quad (157)$$

ACKNOWLEDGMENTS

S. G., F. C., and R. B. acknowledge support funding of the French National Research Agency ANR through project 'Dy-prosome' (ANR-21-CE45-0032-02).

[1] H. Risken, The Fokker-Planck Equation (Springer-Verlag, Berlin, 1989).

[2] C. Gardiner, Stochastic Methods: A Handbook for the Natural and Social Sciences (Springer Verlag, Berlin 2009).

[3] W. T. Coffey, Yu. P. Kalmykov, and J. P. Waldron, The Langevin Equation (World Scientific, Singapore, 2004).

[4] A. T. Silva, E. K. Lenzi, L. R. Evangelista, M. K. Lenzi, H. V. Ribeiro, and A. A. Tateishi, J. Math. Phys. 52, 083301 (2011).

[5] A. G. Cherstvy, R. Metzler, Phys. Chem. Chem. Phys. 15, 20220 (2013).

[6] R. Metzler, J.-H. Jeon, A. G. Cherstvy and E. Barkai, Phys. Chem. Chem. Phys. 16, 24128 (2014).

[7] A. G. Cherstvy, D. Vinod, E. Aghion, A. V. Chechkin, and R. Metzler, New J. Phys. 19, 063045 (2017).

[8] W. Wang, A. G. Cherstvy, X. Liu, and R. Metzler, Phys. Rev. E 102, 012146 (2020).

[9] S. Ritschel, A. G. Cherstvy, and R. Metzler, J. Phys. Complex. 2, 045003 (2021).

[10] D. Vinod, A. G. Cherstvy, W. Wang, R. Metzler, and I. M. Sokolov, Phys. Rev. E 105, L012106 (2022).

[11] D. Vinod, A. G. Cherstvy, R. Metzler, and I. M. Sokolov, Phys. Rev. E 106, 034137 (2022).

[12] H. V. Ribeiro, A. A. Tateishi, E. K. Lenzi, R. L. Magin, and M. Perc, Communications Physics 6, 244 (2023).

[13] N. Leibovich and E. Barkai, Phys. Rev. E 99, 042138 (2019).

[14] E. Aghion, D. A. Kessler, and E. Barkai, Phys. Rev. Lett. 122, 010601 (2019).

[15] E. Aghion, D. A. Kessler, and E. Barkai, Chaos, Solitons and Fractals 138, 109890 (2020).

[16] S. Giordano, F. Cleri, and R. Blossey, Phys. Rev. E 107, 044111 (2023).

[17] A. N. Drozdov, Physica A 196, 258-282 (1993).

[18] A. N. Drozdov, Physica A 196, 283-312 (1993).

[19] A. N. Drozdov, Phys. Rev. Lett. 75, 4342-4345 (1995).

[20] G. H. Weiss and M. Gitterman, Phys. Rev. E 51, 122-125 (1995).

[21] A. N. Drozdov and M. Morillo, Phys. Rev. Lett. 77, 3280-3283 (1996).

[22] M. Donoso, J. J. Salgado and M. Soler, J. Phys. A: Math. Gen. 32, 3681-3695 (1999).

[23] R. Constantinescu, N. Costanzino, A. L. Mazzucato, and V. Nistor, J. Math. Phys. 51, 103502 (2010).

[24] W. Cheng, N. Costanzino, J. Liechty, A. Mazzucato, and V. Nistor, SIAM J. Financial Math. 2, 901-934 (2011).

[25] A. Bilal, J. Math. Phys. 61, 061517 (2020).

[26] K. V. Kumar, S. Ramaswamy, and M. Rao, Active elastic dimers: Self-propulsion and current reversal on a featureless track, Phys. Rev. E 77, 020102(R) (2008).

[27] A. Baule, K.V. Kumar, S. Ramaswamy, J. Stat. Mech. 2008 P11008 (2008).

[28] R. Blossey, H. Schiessel, J. Phys. A: Math. Theor. 52, 085601 (2019).

[29] J. Bauermann, B. Lindner, BioSystems 178, 25-31 (2019).

[30] Z. Zhu, G. Ren, X. Zhang, J. Ma, Chaos, Solitons and Fractals 151, 111203 (2021).

[31] Q. Liu, Y. Jia, Phys. Rev. E 70, 041907 (2004).

[32] D. Frigola, L. Casanellas, J. M. Sancho, M. Ibañes, PLoS ONE 7, e31407 (2012).

[33] F. Manca, P.-M. Déjardin, and S. Giordano, Ann. Phys. (Berlin) 528, 381 (2016).

[34] S. Giordano, Eur. Phys. J. B 92, 174 (2019).

[35] S. Giordano, Phys. Rev. E 103, 052116 (2021).

[36] G. T. Landi and M. J. de Oliveira, Phys. Rev. E 89, 022105 (2014).

[37] P. L. Palla, G. Patera, F. Cleri, and S. Giordano, Physica Scripta 95, 075703 (2020).

[38] A. Fuchs, C. Herbert, J. Rolland, M. Wächter, F. Bouchet, J. Peinke, Phys. Rev. Lett. 129, 034502 (2022).

[39] B. Birnir, The Kolmogorov-Obukhov theory of turbulence (Springer, New York, 2013).

[40] W. Coffey, Adv. Chem. Phys. 63, 69 (1985).

[41] B. Øksendal, Stochastic differential equations (Springer Berlin, 2003).

[42] K. Itô, Nagoya Math. J. 1, 35 (1950).

[43] D. L. Fisk, Quasi-martingales and stochastic integrals (Kent: Research monograph, Kent State University, 1963).

[44] R. L. Stratonovich, SIAM J. Control Optim. 4, 362 (1966).

[45] P. Hänggi, H. Thomas, Phys. Rep. 88, 207 (1982).

[46] Yu. L. Klimontovich, Statistical Theory of Open Systems (Kluver Academic, Dordrecht, 1995).

[47] I. M. Sokolov, Chem. Phys. 375, 359 (2010).

[48] S. I. Denisov, W. Horsthemke, P. Hänggi, Eur. Phys. J. B 68, 567-575 (2009).

[49] S. I. Denisov, A.N. Vitrenko, W. Horsthemke, Phys. Rev. E 68, 046132 (2003).

[50] V. Méndez, S.I. Denisov, D. Campos, W. Horsthemke, Phys.
$$D_1 = 2 + \frac{1}{10}x - \frac{17}{50}x^2 - \frac{137}{1500}x^3 - \frac{307}{30000}x^4 - \frac{91}{180000}x^5 + O(x^6), \quad (154)$$

$$D_2 = \frac{3}{5} + \frac{2}{75}x - \frac{10007}{15000}x^2 - \frac{1621}{7500}x^3 + \frac{24527}{875000}x^4 + \frac{1524793}{52500000}x^5 + O(x^6), \quad (155)$$

$$D_3 = -\frac{2599}{1875} - \frac{1407}{3125}x + \frac{57521}{328125}x^2 + \frac{1992491}{26250000}x^3 + \frac{9657743}{78750000}x^4 + \frac{795208357}{11812500000}x^5 + O(x^6), \quad (156)$$

$$D_4 = -\frac{165647}{656250} - \frac{7079}{93750}x + \frac{45878477}{39375000}x^2 + \frac{347956369}{590625000}x^3 + \frac{741347219}{47250000000}x^4 - \frac{45979054151}{8859375000000}x^5 + O(x^6), \quad (157)$$

ACKNOWLEDGMENTS

S. G., F. C., and R. B. acknowledge support funding of the French National Research Agency ANR through project 'Dy-prosome' (ANR-21-CE45-0032-02).

[1] H. Risken, The Fokker-Planck Equation (Springer-Verlag, Berlin, 1989).

[2] C. Gardiner, Stochastic Methods: A Handbook for the Natural and Social Sciences (Springer Verlag, Berlin 2009).

[3] W. T. Coffey, Yu. P. Kalmykov, and J. P. Waldron, The Langevin Equation (World Scientific, Singapore, 2004).

[4] A. T. Silva, E. K. Lenzi, L. R. Evangelista, M. K. Lenzi, H. V. Ribeiro, and A. A. Tateishi, J. Math. Phys. 52, 083301 (2011).

[5] A. G. Cherstvy, R. Metzler, Phys. Chem. Chem. Phys. 15, 20220 (2013).

[6] R. Metzler, J.-H. Jeon, A. G. Cherstvy and E. Barkai, Phys. Chem. Chem. Phys. 16, 24128 (2014).

[7] A. G. Cherstvy, D. Vinod, E. Aghion, A. V. Chechkin, and R. Metzler, New J. Phys. 19, 063045 (2017).

[8] W. Wang, A. G. Cherstvy, X. Liu, and R. Metzler, Phys. Rev. E 102, 012146 (2020).

[9] S. Ritschel, A. G. Cherstvy, and R. Metzler, J. Phys. Complex. 2, 045003 (2021).

[10] D. Vinod, A. G. Cherstvy, W. Wang, R. Metzler, and I. M. Sokolov, Phys. Rev. E 105, L012106 (2022).

[11] D. Vinod, A. G. Cherstvy, R. Metzler, and I. M. Sokolov, Phys. Rev. E 106, 034137 (2022).

[12] H. V. Ribeiro, A. A. Tateishi, E. K. Lenzi, R. L. Magin, and M. Perc, Communications Physics 6, 244 (2023).

[13] N. Leibovich and E. Barkai, Phys. Rev. E 99, 042138 (2019).

[14] E. Aghion, D. A. Kessler, and E. Barkai, Phys. Rev. Lett. 122, 010601 (2019).

[15] E. Aghion, D. A. Kessler, and E. Barkai, Chaos, Solitons and Fractals 138, 109890 (2020).

[16] S. Giordano, F. Cleri, and R. Blossey, Phys. Rev. E 107, 044111 (2023).

[17] A. N. Drozdov, Physica A 196, 258-282 (1993).

[18] A. N. Drozdov, Physica A 196, 283-312 (1993).

[19] A. N. Drozdov, Phys. Rev. Lett. 75, 4342-4345 (1995).

[20] G. H. Weiss and M. Gitterman, Phys. Rev. E 51, 122-125 (1995).

[21] A. N. Drozdov and M. Morillo, Phys. Rev. Lett. 77, 3280-3283 (1996).

[22] M. Donoso, J. J. Salgado and M. Soler, J. Phys. A: Math. Gen. 32, 3681-3695 (1999).

[23] R. Constantinescu, N. Costanzino, A. L. Mazzucato, and V. Nistor, J. Math. Phys. 51, 103502 (2010).

[24] W. Cheng, N. Costanzino, J. Liechty, A. Mazzucato, and V. Nistor, SIAM J. Financial Math. 2, 901-934 (2011).

[25] A. Bilal, J. Math. Phys. 61, 061517 (2020).

[26] K. V. Kumar, S. Ramaswamy, and M. Rao, Active elastic dimers: Self-propulsion and current reversal on a featureless track, Phys. Rev. E 77, 020102(R) (2008).

[27] A. Baule, K.V. Kumar, S. Ramaswamy, J. Stat. Mech. 2008 P11008 (2008).

[28] R. Blossey, H. Schiessel, J. Phys. A: Math. Theor. 52, 085601 (2019).

[29] J. Bauermann, B. Lindner, BioSystems 178, 25-31 (2019).

[30] Z. Zhu, G. Ren, X. Zhang, J. Ma, Chaos, Solitons and Fractals 151, 111203 (2021).

[31] Q. Liu, Y. Jia, Phys. Rev. E 70, 041907 (2004).

[32] D. Frigola, L. Casanellas, J. M. Sancho, M. Ibañes, PLoS ONE 7, e31407 (2012).

[33] F. Manca, P.-M. Déjardin, and S. Giordano, Ann. Phys. (Berlin) 528, 381 (2016).

[34] S. Giordano, Eur. Phys. J. B 92, 174 (2019).

[35] S. Giordano, Phys. Rev. E 103, 052116 (2021).

[36] G. T. Landi and M. J. de Oliveira, Phys. Rev. E 89, 022105 (2014).

[37] P. L. Palla, G. Patera, F. Cleri, and S. Giordano, Physica Scripta 95, 075703 (2020).

[38] A. Fuchs, C. Herbert, J. Rolland, M. Wächter, F. Bouchet, J. Peinke, Phys. Rev. Lett. 129, 034502 (2022).

[39] B. Birnir, The Kolmogorov-Obukhov theory of turbulence (Springer, New York, 2013).

[40] W. Coffey, Adv. Chem. Phys. 63, 69 (1985).

[41] B. Øksendal, Stochastic differential equations (Springer Berlin, 2003).

[42] K. Itô, Nagoya Math. J. 1, 35 (1950).

[43] D. L. Fisk, Quasi-martingales and stochastic integrals (Kent: Research monograph, Kent State University, 1963).

[44] R. L. Stratonovich, SIAM J. Control Optim. 4, 362 (1966).

[45] P. Hänggi, H. Thomas, Phys. Rep. 88, 207 (1982).

[46] Yu. L. Klimontovich, Statistical Theory of Open Systems (Kluver Academic, Dordrecht, 1995).

[47] I. M. Sokolov, Chem. Phys. 375, 359 (2010).

[48] S. I. Denisov, W. Horsthemke, P. Hänggi, Eur. Phys. J. B 68, 567-575 (2009).

[49] S. I. Denisov, A.N. Vitrenko, W. Horsthemke, Phys. Rev. E 68, 046132 (2003).

[50] V. Méndez, S.I. Denisov, D. Campos, W. Horsthemke, Phys.
Rev. E **90**, 012116 (2014).

[51] G. E. Uhlenbeck and L. S. Ornstein, Phys. Rev. **36**, 823 (1930).

[52] M. C. Wang and G. E. Uhlenbeck, Rev. Mod. Phys. **17**, 323 (1945).

[53] R. N. Mantegna, H.E. Stanley, An Introduction to Econophysics (Cambridge University Press, 2000).

[54] J.-P. Bouchaud, M. Potters, Theory of Financial Risk and Derivative Pricing (Cambridge University Press, 2009).

[55] J. C. Hull, Options, Futures, and Other Derivatives (Pearson Education Ltd., 2021).

[56] V. Stojkoski, T. Sandev, L. Basnarkov, L. Kocarev, and R. Metzler, Entropy **22**, 1432 (2020).

[57] T. L. Saaty, Modern nonlinear equations (Dover Publications, New York, 1981).

[58] S. Hapca, J. Crawford, R. Rae, M. Wilson, and I. Young, Biological Control **41**, 223-229 (2007).

[59] S. Hapca, J. W. Crawford, K. MacMillan, M. J. Wilson and I. M. Young, J. Theor. Biol. **248**, 212 (2007).

[60] S. Hapca, J. W. Crawford and I. M. Young, J. R. Soc., Interface **6**, 111 (2009).

[61] S. B. Yuste, E. Abad and K. Lindenberg, Phys. Rev. E **82**, 061123 (2010).

[62] F. Crick, Nature (London) **225**, 420 (1970).

[63] J. Kowall, D. Peak and J. W. Corbett, Phys. Rev. B : Condens. Matter Mater. Phys., 1976, **13**, 477.

[64] A. G. Kesarev and V. V. Kondrat'ev, Phys. Met. Metallogr. **108**, 30 (2009).

[65] A. G. Kesarev, V. V. Kondrat'ev, and I. L. Lomaev, Phys. Met. Metallogr. **116**, 225-234 (2015).

[66] K. S. Fa, Phys. Rev. E **72**, 020101(R) (2005).

[67] M. Kimura, Journal of Applied Probability **1**, 177-232 (1964).

[68] N. S. Goel, S. C. Maitra, and E. W. Montroll, Rev. Mod. Phys. **43**, 231-276 (1971).

[69] C. P. D. Birch, Annals of Botany **83**, 713-723 (1999).

[70] S. Sakanoue, Ecological Modelling **205**, 159-168 (2007).

[71] K. S. Fa, Langevin and Fokker-Planck equations and their generalizations (World Scientific, Singapore, 2018).
Rev. E **90**, 012116 (2014).

[51] G. E. Uhlenbeck and L. S. Ornstein, Phys. Rev. **36**, 823 (1930).

[52] M. C. Wang and G. E. Uhlenbeck, Rev. Mod. Phys. **17**, 323 (1945).

[53] R. N. Mantegna, H.E. Stanley, An Introduction to Econophysics (Cambridge University Press, 2000).

[54] J.-P. Bouchaud, M. Potters, Theory of Financial Risk and Derivative Pricing (Cambridge University Press, 2009).

[55] J. C. Hull, Options, Futures, and Other Derivatives (Pearson Education Ltd., 2021).

[56] V. Stojkoski, T. Sandev, L. Basnarkov, L. Kocarev, and R. Metzler, Entropy **22**, 1432 (2020).

[57] T. L. Saaty, Modern nonlinear equations (Dover Publications, New York, 1981).

[58] S. Hapca, J. Crawford, R. Rae, M. Wilson, and I. Young, Biological Control **41**, 223-229 (2007).

[59] S. Hapca, J. W. Crawford, K. MacMillan, M. J. Wilson and I. M. Young, J. Theor. Biol. **248**, 212 (2007).

[60] S. Hapca, J. W. Crawford and I. M. Young, J. R. Soc., Interface **6**, 111 (2009).

[61] S. B. Yuste, E. Abad and K. Lindenberg, Phys. Rev. E **82**, 061123 (2010).

[62] F. Crick, Nature (London) **225**, 420 (1970).

[63] J. Kowall, D. Peak and J. W. Corbett, Phys. Rev. B : Condens. Matter Mater. Phys., 1976, **13**, 477.

[64] A. G. Kesarev and V. V. Kondrat'ev, Phys. Met. Metallogr. **108**, 30 (2009).

[65] A. G. Kesarev, V. V. Kondrat'ev, and I. L. Lomaev, Phys. Met. Metallogr. **116**, 225-234 (2015).

[66] K. S. Fa, Phys. Rev. E **72**, 020101(R) (2005).

[67] M. Kimura, Journal of Applied Probability **1**, 177-232 (1964).

[68] N. S. Goel, S. C. Maitra, and E. W. Montroll, Rev. Mod. Phys. **43**, 231-276 (1971).

[69] C. P. D. Birch, Annals of Botany **83**, 713-723 (1999).

[70] S. Sakanoue, Ecological Modelling **205**, 159-168 (2007).

[71] K. S. Fa, Langevin and Fokker-Planck equations and their generalizations (World Scientific, Singapore, 2018).