Feynman-Kac equation for anomalous processes with space- and time-dependent forces

TL;DR

This paper derives a generalized Feynman-Kac equation for anomalous diffusive processes influenced by space- and time-dependent forces, extending the mathematical framework to better model complex biological and physical systems.

Contribution

It provides the first derivation of the Feynman-Kac equation for anomalous processes with space- and time-dependent forces using subordination, and extends it to functionals relevant for stochastic thermodynamics.

Findings

Derived the Feynman-Kac equation for complex anomalous processes

Extended the equation to include time-dependent functionals

Obtained exact results on work fluctuations in a non-equilibrium model

Abstract

Functionals of a stochastic process Y(t) model many physical time-extensive observables, e.g. particle positions, local and occupation times or accumulated mechanical work. When Y(t) is a normal diffusive process, their statistics are obtained as the solution of the Feynman-Kac equation. This equation provides the crucial link between the expected values of diffusion processes and the solutions of deterministic second-order partial differential equations. When Y(t) is an anomalous diffusive process, generalizations of the Feynman-Kac equation that incorporate power-law or more general waiting time distributions of the underlying random walk have recently been derived. A general representation of such waiting times is provided in terms of a L\'evy process whose Laplace exponent is related to the memory kernel appearing in the generalized Feynman-Kac equation. The corresponding anomalous processes have been shown to capture nonlinear mean square displacements exhibiting crossovers between different scaling regimes, which have been observed in biological systems like migrating cells or diffusing macromolecules in intracellular environments. However, the case where both space- and time-dependent forces drive the dynamics of the generalized anomalous process has not been solved yet. Here, we present the missing derivation of the Feynman-Kac equation in such general case by using the subordination technique. Furthermore, we discuss its extension to functionals explicitly depending on time, which are relevant for the stochastic thermodynamics of anomalous diffusive systems. Exact results on the work fluctuations of a simple non-equilibrium model are obtained. In this paper we also provide a pedagogical introduction to L\'evy processes, semimartingales and their associated stochastic calculus, which underlie the mathematical formulation of anomalous diffusion as a subordinated process.

Figures

Equations

$$W(t)={\int }_0^{t}U(Y(r),r)\mathrm{d}r,$$

$$\hat{P}(p,y,t)=⟨e^{ipW(t)}\delta (y-Y(t))⟩$$

$$\frac{\partial }{\partial t}\hat{P}(p,y,t)$$

$$\dot{Y}(t)=F(y,t)+\sigma (y,t)\xi (t),$$

$$Y$$

$$Y(t)=y_0+\sum _{j=1}^{N(t)}{\xi }_j.$$

$$X(s)$$

$$Y(t)$$

$$S(t)=\underset{s>0}{\inf }\left\{s:T(s)>t\right\},$$

$$f_Y(y,t)=⟨\delta (y-Y(t))⟩$$

$$={\int }_0^{\infty }\mathrm{d}sh(s,t)f_X(y,s),$$

$$⟨e^{-\lambda T(s)}⟩=e^{-s{\lambda }^{\alpha }}.$$

$$⟨e^{-\lambda T(s)}⟩=e^{-s\Phi (\lambda )},$$

$$\widetilde{\mathrm{MSD}}(\lambda )$$

$$\dot{X}(s)$$

$$G[k(s)]=⟨e^{-{\int }_0^{\infty }k(s)\eta (s)\mathrm{d}s}⟩=e^{-{\int }_0^{\infty }\Phi (k(s))\mathrm{d}s},$$

$$Y=\sum _{j=1}^n{X}_j^{(n)},$$

$${\varphi }_Y(k)$$

$$Y(t)$$

$$\Psi (k,t)=\ln ⟨\exp [ikY(t)]⟩.$$

$$Y(m)$$

$$Y\left(m\right)$$

$$\Psi (k,m)$$

$$\Psi \left(k,\frac{m}{n}\right)=\frac{m}{n}\Psi (k,1).$$

$${\varphi }_Y(k,t)$$

$$G(k)$$

$${\int }_{\mathbb{R}/\{0\}}\mathrm{Max}(|y{|}^2,1)\Pi (\mathrm{d}y)<\infty .$$

$$G(k)$$

$$Y(t)=\sum _{j=1}^{N(t)}{\xi }_j,$$

Full text

I Introduction

In experimental applications one typically measures physical observables $W$, whose time evolution is determined by the underlying dynamics of the system, which is described by some stochastic process $Y$. Such time-extensive quantities are naturally defined as functionals of the process $Y$ in the form:

[TABLE]

where $U(x,t)$ is some prescribed arbitrary function. If $Y$ is a normal diffusive process, these functionals have been employed to model many different physical phenomena by choosing the function $U$ suitably either with an explicit time dependence or without it. For instance, in the linear case $U(x,t)=x$, with $Y$ interpreted as a particle’s velocity, $W$ represents its position and the $Y$–$W$ representation simply describes the stochastic evolution of the system in the phase space risken1989fokker . If instead we choose $U(x,t)=\delta(x)$ and $U(x,t)=\Theta(x)$, $W$ stands for the local and occupation time respectively darling1957occupation ; agmon1984residence ; newman1998diffusive ; berezhkovskii1998residence ; dhar1999residence ; godreche2001statistics ; majumdar2002local ; blanco2003invariance ; mazzolo2004properties ; barkai2006residence ; grebenkov2007residence ; dumonteil2016residence . Other relevant choices are $U(x,t)=x^{2}$ with Eq. (1) interpreted as a one-dimensional spatial integral, in which case $W$ is interpreted as the variance of a fluctuating interface majumdar2005brownian , and $U(x,t)=e^{-\beta\,x}$ with $\beta$ a real positive parameter, which describes the dynamics of integrated stock prices of the Black-Scholes type yor2012exponential . Another important class of functionals has been introduced in the context of the stochastic thermodynamics of driven small scale systems. Specifically, one is interested in the statistical properties of the accumulated mechanical work done by the system when a non-equilibrium driving is imposed by a time-dependence in a potential $V(x,q(t))$ via some prescribed time-dependent protocol $q(t)$. In such a scenario, assuming $Y$ to be the position coordinate, the mechanical work is defined by Eq. (1) with $U(x,t)=\partial_{q}V(x,q(t))\dot{q}(t)$ jarzynski1997nonequilibrium ; sekimoto1997kinetic ; sekimoto1998langevin ; sekimoto2010stochastic ; seifert2012stochastic .

In order to obtain the probability density function (PDF) of $W$ one usually considers quantities of the form:

[TABLE]

for a given initial condition $Y(0)=y_{0}$. We note that $\widehat{P}$ is the Fourier-transform of the joint PDF of the processes $W$ and $Y$, such that, if one could compute it, the marginal PDF of $W$ would be obtained straightforwardly by making its Fourier inverse transform and subsequently integrating it over all $y$. Equivalently to the direct evaluation of the expected value in Eq. (2), $\widehat{P}$ can be obtained by solving the Feynman-Kac (FK) equation majumdar2005brownian :

[TABLE]

where we introduce the most general Fokker-Planck operator for a space-time dependent force and diffusion coefficient: $\mathcal{L}_{\rm FP}(y,t)=-\frac{\partial}{\partial y}F(y,t)+\frac{1}{2}\frac{\partial^{2}}{\partial y^{2}}\sigma^{2}(y,t)$. In the linear case $U(x)=x$, Eq. (3) maps directly onto the Klein-Kramers equation for the joint position-velocity PDF of a Brownian particle risken1989fokker . The relevance of the FK Equation (3) is motivated by the fact that it allows the calculation of expected values over stochastic trajectories $Y$ of the type of Eq. (2) in terms of solutions of second-order partial differential equations, i.e., of non-random equations, and vice versa. In the conventional setting, i.e., that of Eq. (3), the dynamics of $Y(t)$ is described by the Langevin equation:

[TABLE]

where $\xi(t)$ is a Gaussian white noise with zero mean $\left<\xi(t)\right>=0$ and covariance $\left<\xi(t)\xi(t^{\prime})\right>=\delta(t-t^{\prime})$ and the Itô-convention is assumed for the multiplicative term $\sigma(y,t)\xi(t)$ (Appendix A). Note that the FK equation contains as a special case the Fokker-Planck equation to which it reduces by setting $p=0$ in Eqs. (2, 3). Thus, Eq. (3) is the key method to derive the full statistics of a wide range of phenomena modeled by the diffusive dynamics of Eq. (4) risken1989fokker .

In recent years, an intense effort has been dedicated to derive generalizations of both the FK and Klein-Kramers equation, that extend beyond the normal diffusive regime into the anomalous one. First results were obtained either by substituting $\xi(t)$ with a Lévy noise, such that $Y$ describes Lévy flight type dynamics, peseckis1987statistical ; fogedby1994levy ; jespersen1999levy ; lutz2001fractional ; eliazar2003levy or by directly introducing temporal memory integral terms manifest in time fractional operators into the ordinary FK and Klein-Kramers equations, thus accounting for the non-Markovian effects often characterizing anomalous diffusive processes on a purely phenomenological level metzler2000generalized1 ; metzler2000subdiffusive ; metzler2000generalized2 ; barkai2000fractional ; metzler2002superdiffusive ; zoia2012discrete ; fa2013generalized . These fractional FK equations have been successfully used to model, e.g., the dynamics of migrating epithelial cells dieterich2008anomalous and the advection of a fluid particle in turbulence baule2006investigation . However, the relation between such equations and the underlying stochastic dynamics is often not clear eule2007langevin ; eule2012langevin . Thus, more systematic approaches have been adopted, which explicitly assume the process $Y$ to represent a continuous time random walk (CTRW) with jump lengths and waiting times drawn from independent distributions montroll1965random ; metzler2000random . Specifically, in Friedrich2006Anomalous ; Friedrich2006Exact , starting from a random walk description of the CTRW in phase space, a Klein-Kramers equation containing a fractional substantial derivative, which generalizes the ordinary material derivative through the inclusion of explicit retardation effects, was derived. In turgeman2009fractional a fractional FK equation with the same fractional derivative was derived within a similar random walk description of CTRWs with power-law distributed waiting times. Extensions of this approach to space- and space-time-dependent forces have also been discussed in carmi2010distributions ; carmi2011fractional , as well as to inhomogeneous media in shkilev2012equations ; shkilev2016feynman .

Even if these equations are systematic extensions of the conventional FK formula, they do not establish the same correspondence between some anomalous stochastic dynamics and the solutions of fractional partial differential equations as in the conventional picture of Eqs. (3, 4). Instead of using a Langevin-type equation to describe the dynamics of the underlying CTRW, the time evolution of $\widehat{P}$ is derived directly by means of a generalized master equation. Such full correspondence has been established only recently in our work cairoli2015anomalous by using a general representation of the CTRW in terms of a random time change (also called subordination) of a normal diffusive process. This approach allows in particular to capture straightforwardly different waiting time distributions of the CTRW by a monotonically increasing Lévy process in an auxiliary time variable. The characteristic Laplace exponent $\Phi$ of the Lévy process is naturally related to the memory kernel appearing in the generalized FK equation. Our FK formula has been recently confirmed in wu2016tempered , within a random walk approach, for the special case of tempered Lévy-stable distributed waiting times.

As shown in cairoli2015anomalous , by employing a variable parametric form of $\Phi$, one can fit the resulting anomalous process to mean-square displacement (MSD) data displaying a nonlinear crossover between, e.g., subdiffusive and normal diffusive scaling regimes. In addition, the quantitative form of the higher-order correlation functions of both the CTRW and its observables are fully specified for general $\Phi$, such that they can be readily compared with the experimental data to assess the nature of the underlying stochastic process. Evidence of such crossover scaling behavior has been found in a large number of recent experiments of diffusion in biophysical systems ranging from migrating and foraging cells dieterich2008anomalous ; selmeczi2005cell ; selmeczi2008cell ; campos2010persistent ; harris2012generalized to macromolecules and living organelles, e.g., mitochondria, inside the cytoplasm caspi2000enhanced ; levi2005chromatin ; brangwynne2007force ; bronstein2009transient ; bruno2009transition ; senning2010actin ; Jeon2011InVivo ; jeon2012anomalous ; weber2012nonthermal ; von2013anomalous ; tabei2013intracellular ; javer2014persistent . Thus, our framework can be applied to a large variety of different systems exhibiting such anomalous diffusive behavior.

Our main purpose in the present manuscript is to extend the derivation of the generalized FK equation for CTRWs with an arbitrary waiting time distribution cairoli2015anomalous to both space- and time- dependent forces and explicit time-dependent functionals. These latter ones in particular, to our knowledge, have not yet been considered in previous works on anomalous diffusion processes and their observables, despite their great importance in the context of the stochastic thermodynamics of small scale systems. Even though a FK equation for space-time-dependent forces acting on the CTRW has already been presented in carmi2011fractional within a random walk approach for power-law distributed waiting times, its extension to arbitrary distributions as expressed in the subordination approach of cairoli2015anomalous has not been presented so far. Thus, we here provide the missing link to a comprehensive understanding of functionals of anomalous processes. Specifically, our proposed equations will allow to model the effect of non-equilibrium work protocols on biophysical systems exhibiting complex anomalous diffusion. An additional purpose is to provide a largely self-contained and pedagogical introduction to Lévy processes, semimartingales and their stochastic calculus, which is necessary to understand the mathematical framework underlying the description of anomalous processes in terms of subordination.

The remainder of this paper is organized as follows. In Sec. II we review the definition of the CTRW model with arbitrarily distributed waiting times and its representation in the diffusive limit as a subordinated stochastic process, whose dynamics is described by coupled Langevin equations. In Sec. III we provide the mathematical fundamentals that are necessary to manipulate this representation formally. The stochastic calculus of such subordinated processes is briefly discussed. In Sec. IV we use the appropriate form of Ito’s lemma to derive generalized FK equations. For pedagogical reasons, we first treat the case of a space-dependent force and time-independent functional as in cairoli2015anomalous . We then extend the method to space- and time-dependent forces and time-dependent functionals. In Sec. V we apply our results to study the accumulated mechanical work fluctuations in a simple non-equilibrium model with anomalous dynamics. Finally, in Sec. VI we provide some final remarks on open questions and future work.

II Anomalous processes with general waiting times

The discussion of random walks with arbitrarily distributed waiting times goes back to the seminal work by Montroll and Weiss montroll1965random . In this picture, the random walk is defined as a renewal process, where the walker selects the jump lengths as identically and independently distributed (i.i.d.) random variables (RVs). The waiting times between each jumps are also i.i.d. RVs with a distribution that is possibly correlated with the jump length one. In this paper we generally assume that waiting times and jump lengths are uncorrelated. For a discussion of the correlated case, we refer to tejedor2010anomalous ; magdziarz2012correlated ; magdziarz2012langevin ; schulz2013correlated ; de2013flow ; liu2013continuous ; magdziarz2013asymptotic . A natural parametrisation of such a random walk is obtained in terms of the number $n$ of jumps performed. If we call $\xi_{j}$ the amplitude of the jump occurring at the $j$th step and by $\eta_{j}$ the waiting time between the $(j-1)$th and $j$th jumps, the position $Y$ and the elapsed time $T$ are given by summing all such $n$ RVs:

[TABLE]

where $y_{0}$ denotes the initial position. Rather than a parametrisation in terms of the discrete variable $n$, it is usually preferable to describe the position coordinate in terms of a continuous time variable $t$. In Eq. (5) we see that $T$ and $n$ are complementary variables, i.e., either one considers $n$ to be a fixed (integer) number, in which case $T$ is a RV, or one considers directly $n$ to be a RV that gives the number of jumps in a time interval $[0,t]$, where $t$ is the elapsed physical time. If one adopts the latter viewpoint, $n$ becomes the stochastic process $N(t)$, which is defined formally as $N(t)=\max{\{n\geq 0:T(n)\leq t\}}$, and the position variable can be written as below:

[TABLE]

Now the waiting time statistics are contained in $N(t)$. Consequently, there are two main methods to obtain the statistics of $Y(t)$ from this Montroll-Weiss random walk picture: (i) One can formulate a generalized master equation for the PDF $f_{Y}(y,t)\!=\!\langle\delta(y-Y(t))\rangle$ directly from Eq. (6). The master equation is then further approximated on a diffusive time and spatial scale leading to Fokker-Planck equations with fractional time derivatives, which describe the time evolution of $f_{Y}(y,t)$ metzler1998fractional ; metzler1998anomalous ; metzler1999deriving ; metzler1999anomalousPRL ; metzler1999transport ; metzler2000random ; barkai2000continuous . (ii) The diffusive limit can already be considered on the level of Eq. (5), thus leading to a coupled set of Langevin equations describing the stochastic process $Y(t)$ fogedby1994langevin ; baule2005joint ; weron2008modeling . The resulting Fokker-Planck equation for $f_{Y}(y,t)$ is equivalent to that obtained by using approach (i) baule2005joint ; magdziarz2007fractional ; magdziarz2008equivalence ; henry2010fractional .

In the following, we focus on (ii) and provide a pedagogical introduction to the mathematical framework that is needed to describe anomalous diffusive systems within this approach. The key step is to take a continuum limit in the number of steps: $N(t)\to S(t)$ meerschaert2011fractional . In such a continuum limit Eqs. (5) become:

[TABLE]

where $s$ is now interpreted as an auxiliary or operational time variable. The position coordinate Eq. (6) becomes:

[TABLE]

Thus, one must distinguish the two processes $Y$ and $X$, which are parametrised by the physical and the auxiliary time respectively. The complementary relationship between them is naturally expressed by

[TABLE]

i.e., $S$ is defined as a collection of first passage times. Indeed, this definition ensures that $S$ accounts exactly for the number of steps, such that the total elapsed time, i.e., the sum of the waiting time increments for each of those steps, is equal to $t$. We see that Eq. (9) defines $S$ formally as the inverse process of $T$. Thus, the CTRW $Y(t)$ is naturally defined as $Y(t)=X(S(t))$, i.e., as a time-changed or subordinated process. The mathematical details underlying the subordination concept are addressed in Sec. III.3. An illustrative picture of this representation, compared with the ordinary random walk, i.e., a normal diffusion, is presented in Fig. 1. We note that under the assumption of uncorrelated jump lengths and waiting times the PDF of $Y(t)$ can be expressed in terms of the integral transform:

[TABLE]

where $f_{X}$ is the PDF associated with the process $X(s)$ and $h(s,t)$ that of $S(t)$.

Regarding the waiting time process, a widely studied case is that of $T$ being a one-sided Lévy-stable process of order $0<\alpha\leq 1$, which corresponds to a distribution of the waiting times with power-law tails and diverging first moment. In this specific case, $T$ has the characteristic function:

[TABLE]

Thus, the time-change $S$ in Eq. (9) is an inverse Lévy-stable subordinator with a PDF defined in Laplace space as $\widetilde{h}(s,\lambda)=\int_{0}^{\infty}\mathrm{d}t\,e^{-\lambda t}h(s,t)=\lambda^{\alpha-1}\,e^{-s\,\lambda^{\alpha}}$ baule2005joint . When $X$ describes a pure diffusion process with noise strength $\sigma$, the MSD of $Y$ exhibits a power-law scaling of the same exponent $\alpha$, i.e., ${\rm MSD}(t)=\left<[Y(t)-y_{0}]^{2}\right>=\sigma\,t^{\alpha}/\Gamma(1+\alpha)$. This particular scaling regime is characteristic of a pure subdiffusive system metzler2000random .

However, in realistic situations the MSD does not always exhibit a single power-law scaling. In fact, diffusive systems, whose MSD exhibits possibly multiple crossovers between different scaling regimes, are widely observed selmeczi2005cell ; selmeczi2008cell ; dieterich2008anomalous ; harris2012generalized ; campos2010persistent ; caspi2000enhanced ; levi2005chromatin ; brangwynne2007force ; bronstein2009transient ; bruno2009transition ; senning2010actin ; jeon2012anomalous ; weber2012nonthermal ; von2013anomalous ; tabei2013intracellular ; javer2014persistent . The generalization of the power-law case to account for such more general MSDs is obtained mathematically by choosing $T$ to be a general one-sided strictly increasing Lévy process. Indeed, such a process satisfies the minimal assumptions needed to assure independent and stationary waiting times and causality of $T$. Thus, we specify $T$ by means of its characteristic function, which is given by cont1975financial ; applebaum2009levy :

[TABLE]

with the Laplace exponent $\Phi$ characterizing the jump structure of the waiting times. Different functional forms of $\Phi$ correspond to different distribution laws of the waiting times and of the renewal process $T$. By choosing $\Phi$ suitably, several different waiting time statistics can be captured, i.e., the anomalous process $Y(t)$ can be modeled according to the observed experimental dynamics. If we choose a power law $\Phi(\lambda)=\lambda^{\alpha}$, we recover Eq. (11), i.e., the CTRW case. If instead $\Phi(\lambda)=\lambda$, $T$ is a deterministic drift, $T=s$, and $Y(t)$ reduces to a normal diffusion (Brownian limit) with exponentially distributed waiting times metzler2000random . Details on the mathematical properties of $\Phi$ are discussed in Sec. III.2.

For $X$ given as a normal diffusion, the MSD of the anomalous process $Y$ with general waiting times can be computed straightforwardly by employing Eq. (II). Indeed, the inverse of the process $T$ has the PDF in Laplace space: $\widetilde{h}(s,\lambda)=[\Phi(\lambda)/\lambda]\,e^{-s\,\Phi(\lambda)}$ magdziarz2009langevin . As the first and second moment of $X$ are given by $\langle X(s)\rangle=y_{0}$ and $\langle[X(s)]^{2}\rangle=y^{2}_{0}+\sigma s$ respectively, those of the time-changed process $Y$ read as $\langle\widetilde{Y}(\lambda)\rangle=\int_{0}^{\infty}\widetilde{h}(s,\lambda)\,\langle X(s)\rangle=y_{0}/\lambda$ and $\langle\widetilde{Y}^{2}(\lambda)\rangle=\int_{0}^{\infty}\widetilde{h}(s,\lambda)\,\langle X^{2}(s)\rangle=[y_{0}^{2}+\sigma/\Phi(\lambda)]/\lambda$. Putting these results together, we obtain the MSD in Laplace space:

[TABLE]

We note that for $\Phi(\lambda)=\lambda^{\alpha}$ we recover the single power-law scaling previously discussed. For different choices of $\Phi$ Eq. (13) is able to capture many different scaling behaviors of ${\rm MSD}(t)$ cairoli2015anomalous . For instance, in the case of $T$ given as a tempered Lévy-stable process, i.e., $\Phi(\lambda)=(\mu+\lambda)^{\alpha}-\mu^{\alpha}$ with $\mu\in\mathbb{R}^{+}$, the MSD displays a crossover between a power-law (of exponent $\alpha$) and a normal linear scaling. In the case of a sum of two independent stable distributions with exponents $\alpha_{1}<\alpha_{2}$, i.e., $\Phi(\lambda)=C_{1}\,\lambda^{\alpha_{1}}+C_{2}\,\lambda^{\alpha_{2}}$ with $C_{1},C_{2}\in\mathbb{R}^{+}$, the crossover is between two subdiffusive regimes with power-law scaling of exponents $\alpha_{1}/\alpha_{2}$ respectively for long/small times sandev2015distributed .

Equations (7) allow us to express $X$ and $T$ in terms of Langevin equations, as first formulated by Fogedby fogedby1994langevin :

[TABLE]

with initial conditions $X(0)=Y(0)=y_{0}$ and $T(0)=0$. Even though these steps seem somewhat superfluous, there are various advantages by expressing both $X$ and $T$ in this way. In particular, it allows to easily incorporate forces acting on the random walker during the instantaneous jumps. The case where they affect the dynamics of the walker also during the waiting times is discussed in fedotov2014sub ; cairoli2015langevin . The factorization of the expected values leading to Eq. (II) still holds in the presence of forces that depend only on the position or on the auxiliary time variable, i.e., when Eq. (14)(left) is substituted by $\dot{X}(s)=F(X(s),s)+\xi(s)$. However, in realistic scenarios the force should vary in physical time rather than in the auxiliary one, which is only a formal construct to simplify the mathematical description heinsalu2007use . Instead, the correct Langevin equation substituting Eq. (14)(left) is given as $\dot{X}(s)=F(X(s),T(s))+\xi(s)$ heinsalu2007use ; magdziarz2008equivalence ; weron2008modeling ; magdziarz2009stochastic ; magdziarz2009langevin ; eule2009subordinated ; heinsalu2009fractional ; henry2010fractional . Now the factorization in Eq. (II) breaks down. Indeed, the evolution of the position in the auxiliary time depends on the statistics of both the jumps lengths and waiting times. Therefore, introducing space- and (physical) time-dependent forces naturally induces a coupling between the two corresponding processes.

By specifying the properties of $X$ and $T$, we can capture the fluctuation properties of a large variety of physical systems. In the following, we generally assume $\xi(s)$ to represent a white Gaussian noise with $\left\langle\xi(s)\right\rangle=0$ and $\left\langle\xi(s)\xi(s^{\prime})\right\rangle=\delta(s-s^{\prime})$, such that $X$ is a normal diffusive process in the auxiliary time $s$. For generality, we also consider a multiplicative noise strength $\sigma(x)$, which can capture, e.g., the effects of geometrical confinement (see lau2007state and references therein). Thus, we consider the following set of coupled Langevin equations:

[TABLE]

where the functions $F(x,t)$ and $\sigma(x)$ satisfy standard conditions revuz1999continuous , we adopt the Itô prescription for the multiplicative term (Appendix A) and we assume the same initial conditions of Eqs. (14). The stochastic process $Y$ describing the dynamics of the position coordinate is again obtained by subordination, i.e., by Eqs. (8, 9).

We note that instead of defining $T$ directly by Eq. (12), it is convenient to specify the underlying noise process $\eta(s)$ in Eq. (15)(right) by means of its characteristic functional van1992stochastic ; cont1975financial :

[TABLE]

which specifies the statistics of the whole noise trajectory. Recalling Eq. (7), we recover the characteristic function $\left<e^{-\lambda T(s)}\right>$ from Eq. (16) by setting $k(s^{\prime})=\lambda\,\Theta(s-s^{\prime})$. This result, together with Eq. (15)(right), elucidates the renewal nature of the process $T$. Such a process is indeed expressed as a sum over waiting time increments $\Delta t=\int_{0}^{\Delta s}\eta(\tau)\mathop{}\!\mathrm{d}\tau$ over a small time step $\Delta s$ of characteristic function $\left<e^{-\lambda\Delta t}\right>=e^{-\Delta s\Phi(\lambda)}$, which can be used to simulate the process $Y(t)$ within a suitable discretisation scheme kleinhans2007continuous . Specifying the noise $\eta$ by means of its characteristic functional Eq. (16) also renders the full multi-point statistics of $T$ easily accessible. In the next section we briefly review some fundamental mathematical notions regarding Lévy processes, semimartingales and time-changed processes. In particular, we characterize the function $\Phi$, which specifies the waiting time statistics of the underlying random walk.

III Mathematical basics: Lévy processes, subordinators and time-changed processes

In this section, we provide a comprehensive and accessible review of the theory of Lévy processes, subordinators and time-changed processes, that are necessary to formulate and work with Langevin equations of anomalous diffusive processes. For more details beyond our present discussion we refer to the monographs cont1975financial ; applebaum2009levy .

III.1 Lévy processes

A stochastic process $Y(t)$ for $t\geq 0$ and initial condition $Y(0)=y_{0}$ is a Lévy process if the conditions below hold:

$Y(0)=y_{0}=0$ almost surely (a.s.), i.e., for each of its different realizations. 2. 2.

$Y(t)$ has independent increments, i.e., $\forall\,n\geq 2$ and for each partition $0\leq t_{0}<t_{1}<\ldots<t_{n}\leq t$ the RVs $\{Y(t_{j})-Y(t_{j-1})\}_{j=1,\ldots,n}$ are independent. 3. 3.

$Y(t)$ has stationary increments, meaning that for all $0\leq t_{1}<t_{2}\leq t$ the RV $Y(t_{2})-Y(t_{1})$ has the same distribution as $Y(t_{2}-t_{1})$. Note that, if 1 is not satisfied, it would instead depend on $Y(t_{2}-t_{1})-y_{0}$. 4. 4.

The trajectories of $Y(t)$ are càdlàg, i.e., right continuous with left limits.

If one restricts the conditions 2, 4 by assuming Gaussian distributed increments and continuous trajectories respectively, one recovers ordinary Brownian motion. Moreover, as a consequence of (ii), $Y$ is infinitely divisible $\forall\,t\geq 0$.

The notion of infinitely divisibility characterizes a RV $Y$ that can be expressed as a sum of different i.i.d. RVs $\{X_{j}\}$. Specifically, we assume $Y$ to have PDF $f_{Y}$ and corresponding characteristic function: $\phi_{Y}(k)=\left\langle e^{i\,k\,Y}\right\rangle=\int_{-\infty}^{+\infty}e^{i\,k\,x}\,f_{Y}(x)\mathop{}\!\mathrm{d}x$. If $\forall\,n\in\mathbb{N}$ there exist i.i.d RVs $X^{(n)}_{1},\ldots,X^{(n)}_{n}$ with PDF $f_{X}$ and characteristic function $\phi_{X}$ (uniquely defined), such that the following relation holds (in distribution):

[TABLE]

then $Y$ is said to be infinitely divisible. Thus, their characteristic functions are related by the following equation:

[TABLE]

Note that the factorization of the average is due to the independence of the RVs $X_{j}^{(n)}$ and that Eq. (18) represents a necessary and sufficient condition for $Y$ to be infinitely divisible applebaum2009levy , i.e., it can be used as a criterion to asses the infinitely divisibility of a given RV.

Considering now the Lévy process $Y$ in continuous time, we find that $\forall\,n\in\mathbb{N}$ and $\forall\,t\geq 0$ it can be rewritten as:

[TABLE]

where $\{X^{(n)}_{j}\}$ are i.i.d. RVs because of 2-3. Therefore, such property uniquely relates the characteristic function of Y to that of a general infinitely divisible RV. Indeed, let us introduce the auxiliary function:

[TABLE]

By suitably adapting Eq. (19), we can write for arbitrary $m,n\in\mathbb{N}$:

[TABLE]

Further exploiting the property 2, Eqs. (21a, 21b) imply that $Y(m)=m\,Y(1)$ and $Y(m)=n\,Y\!\left(\frac{m}{n}\right)$ (in distribution), such that $\Psi(k,m)$ can be computed exactly. Specifically, we obtain the two equivalent equations:

[TABLE]

Here, the factorization of the ensemble average is allowed by the independence of the increments. Consequently, the equations $\Psi(k,m)=m\,\Psi(k,1)$ and $\Psi(k,m)=n\,\Psi\!\left(k,m/n\right)$ hold simultaneously. If we combine them, we obtain:

[TABLE]

This relation expresses the characteristic function of Y at the finite time $t=m/n$ in terms of its value at time $t=1$. As it is satisfied for every integer $m,n\,$, it also holds for any real positive $t$, i.e., $\Psi\!\left(k,t\right)=t\,\Psi(k,1)$. Moreover, according to Eq. (19), $Y(1)$ is an infinitely divisible RV. Therefore, a Lévy process $Y(t)$ has the characteristic function:

[TABLE]

where we set $G(k)=\Psi(k,1)$, with $\Psi(k,1)$ being fully specified as the logarithm of the characteristic function of a general infinitely divisible RV. Such quantity, along with the characteristic function of a general Lévy process according to Eq. (24), is uniquely characterized by means of the Lévy-Khintchine representation. Introducing parameters $b\in\mathbb{R}$ and $\sigma\geq 0$, the Lévy-Khintchine representation states that the characteristic function of any infinitely divisible RV is of the form $\phi(k)=\exp{[G(k)]}$ with the function $G$ defined below applebaum2009levy :

[TABLE]

with $\bm{1}_{A}(y)=1$ for $y\in A$ or $\bm{1}_{A}(y)=0$ otherwise. In Eq. (25) $\Pi$ is a so called Lévy measure, i.e., a probability measure satisfying the following condition:

[TABLE]

Thus, any Lévy process $Y$ is uniquely characterized by the triplet $(b,\sigma,\Pi)$, which determines its characteristic function trough Eqs. (24, 25). Two examples of Lévy processes are of fundamental importance:

• •

Brownian motion with drift. In this case, the characteristic function is given by Eq. (24) with the specification:

[TABLE]

As we can write $G$ as below:

[TABLE]

where we define $G_{X}(k)=\exp{[i\,k\,b/n-k^{2}\sigma^{2}/(2\,n)]}$, we deduce that $G(k)$ is the characteristic function of an infinitely divisible RV. The RVs $Y^{(n)}_{j}$ in Eq. (17) are Gaussian distributed with mean $b/n$ and variance $\sigma^{2}/n$.

• •

Compound Poisson process. A Compound Poisson process $Y(t)$ on the interval $[0,t]$ is defined as

[TABLE]

where the $\xi_{j}$ are i.i.d. RVs with law $f_{\xi}$ and $N(t)$ is a Poisson process characterized by an intensity $\lambda$, i.e., $p(N(t)=n)=\exp(-\lambda t)[(\lambda t)^{n}/n!]$. This represents a pure jump process where jumps of length $\xi_{j}$, drawn as i.i.d. RVs, occur at time points that are spaced by an exponentially distributed waiting time. Within this picture, $\lambda$ represents the average number of jumps per unit time. Comparing it with Eq. (6), we note that the Compound Poisson process can be regarded as the simplest renewal process. Its characteristic function can be calculated in a straightforward way by conditioning on the number of jumps:

[TABLE]

Therefore, $Y(t)$ is a Lévy process with

[TABLE]

Thus, $G(k)$ is again the characteristic function of an infinitely divisible RV, as Eq. (17) is satisfied by choosing the RVs $Y^{(n)}_{j}$ to have the characteristic function Eq. (31) with $\lambda\to\lambda/n$.

Comparing these two examples with the Lévy-Khintchine Eq. (25) we note that Lévy processes can be interpreted intuitively as consisting of three different contributions: A deterministic drift, a continuous normal diffusion and a discontinuous jump process. The Compound Poisson process represents the simplest example of the jump contribution, where the Lévy measure of the jump amplitudes is just a normalizable PDF: $\Pi(\mathrm{d}y)=\lambda\,f_{\xi}(y)\mathrm{d}y$. The specifications given in Eqs. (25, 26) extend this case to a wider class of jump processes, which have possibly non-normalizable length distribution or may possess an infinite intensity of small jumps gardiner1985handbook . For physical applications, important examples of such more exotic processes are Lévy-stable and tempered Lévy-stable processes samoradnitsky1994stable .

A Lévy-stable process is a Lévy process with stable distributed increments, i.e., $G(k)$ in Eq. (24) is the characteristic function of a stable RV, which constitutes a special case of infinitely divisible RVs. Let us consider a RV $Y$ and $n$ independent of its copies $\{Y_{j}\}_{j=1,\ldots,n}$. If real-valued sequences of parameters $\{c_{n}\}_{n\in\mathbb{N}}$ and $\{d_{n}\}_{n\in\mathbb{N}}$ exist, such that the following relation holds in distribution:

[TABLE]

then $Y$ is called a stable RV. If $d_{n}=0$, then $Y$ is strictly stable. From this definition, it is straightforward to see that (i) $Y$ is infinitely divisible [simply set $X_{j}^{(n)}=(Y_{j}-d_{n}/n)/c_{n}$ in Eq. (17)] and that (ii) the existence of $Y$ represents a generalization of the central limit theorem. Indeed, Eq. (32) equivalently states that the sequences of partial sums $\{S_{n}\}_{n\in\mathbb{N}}$ with $S_{n}=(Y_{1}+\ldots+Y_{n}-d_{n})/c_{n}$ converge in distribution to $Y$. With the choice $c_{n}=\sigma\,\sqrt{n}$ and $d_{n}=n\,m$, this is the ordinary central limit theorem and $Y$ is Gaussian distributed with mean $m$ and variance $\sigma^{2}$. For different choices of $c_{n}$ and $d_{n}$, we obtain instead a generalized central limit theorem gnedenko1954limit . However, the only possible choice to satisfy Eq. (32) is given by $c_{n}=\sigma\,n^{1/\alpha}$, with $0<\alpha\leq 2$, also called index of stability of the stable distribution feller1971introduction . As stable distributions are infinitely divisible, their characteristic function is completely determined by Eq. (25). In particular, we have two possible characteristics: (i) $(b,\sigma^{2},0)$ for $\alpha=2$, implying that $Y$ is Gaussian (mean $b$, variance $\sigma^{2}$) and (ii) $(b,0,\Pi)$ for $\alpha\neq 2$ with $\Pi$ specified by the following formula (for $c_{1},c_{2}\geq 0$ and $c_{1}+c_{2}>0$):

[TABLE]

By suitably changing coordinates in Eq. (25) sato1999levy , we obtain the following characterization of $G(k)$:

[TABLE]

for $\mu\in\mathbb{R}$, $\sigma\geq 0$ and $-1\leq\beta\leq 1$. If $X$ is a symmetric stable RV ($\beta=0$), then the function $G$ reduces to:

[TABLE]

with $\rho=\sigma$ for $0<\alpha<2$ and $\rho=\sigma/\sqrt{2}$ for $\alpha=2$.

Of particular importance for us are one-sided monotonically increasing Lévy processes, which can be used to implement a random time change. These processes are called subordinators.

III.2 Subordinators

We define a subordinator a one-dimensional Lévy process that is a.s. non-decreasing. Thus, if $T(t)$ for $t\geq 0$ is a subordinator, the following properties hold a.s.: (i) $T(t)\geq 0$, $\forall\,t\geq 0$ and (ii) $T(t_{1})\leq T(t_{2})$, $\forall\,t_{1}\leq t_{2}$. According to the discussion in Sec. III.1, its characteristic function is determined by Eqs. (24, 25) for a subclass of characteristic triplets $(b,\sigma,\nu)$ that we need to determine. We first note that, if $X(t)$ is a Brownian motion of variance $\sigma^{2}$, we have: $p(X(t)\geq 0)=1/2=p(X(t)\leq 0)$. For a subordinator instead, we require $p(T(t)<0)=0$ for all times. Thus, a subordinator $T$ cannot have any Gaussian component in its Lévy symbol, i.e., $\sigma=0$ in Eq. (25). In addition, the monotonicity of $T$ implies that no jumps of negative amplitudes nor a negative shift are allowed, thus implying the further conditions: $b\geq 0$ and $\Pi(-\infty,0)=0$. Taking these requirements into account, the characteristic function of a general subordinator $T$ is given as bertoin1999subordinators :

[TABLE]

where one needs to further assume that $\int_{0}^{+\infty}\mathrm{Max}{(y,1)}\,\Pi(\mathop{}\!\mathrm{d}y)<\infty$. $\Phi$ is the Laplace exponent of the subordinator. As suggested in Sec. II, $\Phi$ determines the characteristic functional Eq. (16) of the noise $\eta$ appearing in Eq. (15)(right). We remark that only two parameters define its form, i.e., the characteristics of $T$ are determined by the duplet $(b,\Pi)$. Using Eq. (24) and Jensen’s inequality, one can show that $\Phi(\lambda)$ must be a continuous, non negative, non decreasing and concave function. We also remark that $\Phi(0)=0$. In general, one can prove that $\Phi$ is a Bernstein function schilling2012bernstein ; meerschaert2015relaxation . Specific examples of subordinators are reviewed in the following:

• •

Lévy stable subordinator. A subordinator $T$ is Lévy stable if it has characteristic duplet $(0,\Pi)$ with

[TABLE]

If we substitute it inside Eq. (37), we obtain the following Laplace exponent:

[TABLE]

• •

Tempered Lévy stable subordinator. A subordinator $T$ is tempered Lévy stable if it has characteristic duplet $(0,\Pi)$ with the Lévy measure cont1975financial :

[TABLE]

If we substitute it inside Eq. (37), we obtain the following Laplace exponent:

[TABLE]

where we solved the integrals by employing Eq. (39). The resulting process $T$ interpolates between a Lévy stable (of order parameter $\alpha$) and an exponential process. This can be shown by applying the Tauberian theorems, which relate the long(small)-$T$ limit of its distribution to the $\lambda\to 0$($\infty$) limit of its Laplace transform Eq. (36), or equivalently of the function $\Phi$ derived in Eq. (41). Specifically, for $\lambda\to\infty$ we find: $\Phi(\lambda)\sim\lambda^{\alpha}$, which recovers the case of a Lévy stable subordinator [see Eq. (39)]. For $\lambda\to 0$ instead, we obtain: $\Phi(\lambda)\sim\alpha\,\mu^{\alpha-1}\,\lambda$ , such that we can approximate the characteristic function of $T$ as $\langle e^{-\lambda\,T(t)}\rangle\sim 1-\alpha\,\mu^{\alpha-1}\,t\,\lambda\sim(1+\alpha\,\mu^{\alpha-1}\,t\,\lambda)^{-1}$, which is the Laplace transform of an exponential distribution.

III.3 Semimartingales and the stochastic calculus of time-changed processes

Let us consider a Lévy process $X(t)$ and a subordinator $T(t)$. Thanks to its monotonicity, $T$ can be employed directly as a random parametrisation of time defining a new time-changed process $Y$ trough the relation: $Y(t)=X(T(t))$. Such process is easily shown to still be a Lévy process applebaum2009levy . However, as discussed in Sec. II, this is not the situation arising for CTRWs, where their representation by coupled Langevin equations in the diffusive limit involves the inverse of the process $T$, i.e., the process $S$ defined in Eq. (9). The crucial point is that $S$ is generally not a Lévy process. Rather, it is part of a more general class of processes called semimartingales, which also contain Lévy processes as a special case (see below). An important theorem tracing back to the work of Jacod jacod714calcul states that semimartingales (and thus Lévy processes) subordinated by properly defined time-changes, e.g., the process $S$, are again semimartingales. Thus, when we study anomalous diffusion at the level of the Langevin representation of Eqs. (15), we need to employ the stochastic calculus of semimartingales. Recalling that the trajectories of $S$ are continuous, because $T$ itself is strictly increasing according to Eq. (15)(right) revuz1999continuous ; kobayashi2011stochastic , we can focus only on the subclass of continuous semimartingales.

Let us consider a process $M(t)$ and assume that all the information on $M$ up to a chosen time $s$ is known, i.e., we know $M(s)=m_{s}$. The process $M$ is a martingale if the following relation on its conditional average holds gardiner1985handbook :

[TABLE]

It is instead called a sub-martingale if $\left\langle M(t)|M(s)=m_{s}\right\rangle\geq m_{s}$ or a super-martingale if $\left\langle M(t)|M(s)=m_{s}\right\rangle\leq m_{s}$. For instance, the Brownian motion $B(t)$ is a martingale, as one can easily verify by direct computation of Eq. (42). We define a process $Y$ a semimartingale if the following decomposition holds:

[TABLE]

where $M(t)$ and $A(t)$ are a martingale and a finite variation process with càdlàg paths respectively. We recall that stochastic integration with respect to semimartingales is well defined kunita1997stochastic . For the sake of our discussion, we will only present their Itô formula. In the specific case of $Y$ being a continuous semimartingale, this is given by:

[TABLE]

where $[Y,Y]_{t}$ is the quadratic variation of $Y$ (see Appendix C for a review). The extension of Eq. (44) to a M-dimensional semimartingale $\bm{Z}$ reads as:

[TABLE]

where $[Z^{(i)},Z^{(j)}]_{t}$ is the joint quadratic variation of $Z^{(i)}$, $Z^{(j)}$, which is defined analogously to the quadratic variation by substituting the squared increment in Eq. (126) with the product of the increments of the two processes. We note that both $[Y,Y]_{t}$ and $[Z^{(i)},Z^{(j)}]_{t}$ are continuous increasing processes. The joint one also has finite variation paths kunita1997stochastic . For further properties of semimartingales and their theory of stochastic integration we refer to kunita1997stochastic .

Recalling that both $T(s)$ and $S(t)$ are monotonically non decreasing, we can deduce that $S$ is a process of finite variation (see Appendix B for a justification of this statement). This property, on the one hand, classify $S$ generally as a semimartingale [according to Eq. (43)] and, on the other hand, together with the continuity of its paths, enables us to specify its Itô formula for a general differentiable function $f$ as follows (adapted from Eq. (124) in Appendix B):

[TABLE]

In addition, the monotonicity of the paths of $T$ and $S$ also provides the relation baule2005joint :

[TABLE]

Thus, if we choose $f(S(t))=\Theta(s-S(t))$ in Eq. (46) and we use Eq. (47), we obtain:

[TABLE]

or equivalently in its corresponding differential form cairoli2015anomalous :

[TABLE]

We note that $\dot{S}(t)=\lim_{\Delta t\to 0}[S(t+\Delta t)-S(t)]/\Delta t$ is a shorthand notation to denote an integration with respect to the time-change. With this definition, we can rewrite the coupled Langevin Eqs. (15) as a single time-changed stochastic differential equation kobayashi2011stochastic . To avoid technicalities related to the jumps of $T$, we first neglect the time dependence of the force term, i.e., $F(x,t)\rightarrow F(x)$ in Eq. (15)(left). The general case will be addressed in Sec. IV.2. Thus, we can integrate Eq. (15)(left) directly to obtain ($X(0)=y_{0}$):

[TABLE]

If we now apply directly the time change, we obtain the integrated equation for $Y$:

[TABLE]

In order to proceed, we recall the following two key results valid for time-changed semimartingales. Let $Z$ be a continuous semimartingale and $S$ be given by Eq. (9) for a subordinator $T$. One can prove the following kobayashi2011stochastic :

[TABLE]

where $H$ is any function that can be integrated with respect to $Z$. Applying Eq. (52a) to Eq. (51) yields:

[TABLE]

which can finally be written as a Langevin equation by taking its time derivative:

[TABLE]

This equation directly expresses the evolution of the increments of $Y$ in terms of those of the time-change $S$. The term $\xi(S(t))\dot{S}(t)$ denotes an increment over the time-changed Brownian motion: $\xi(S(t))\dot{S}(t)=\lim_{\Delta t\to 0}[B(S(t+\Delta t))-B(S(t))]/\Delta t$. To justify this result, we recall that increments of the Brownian motion can be written in terms of the noise $\xi$ trough the integral relation: $B(t+\Delta t)-B(t)=\int_{t}^{t+\Delta t}\xi(\tau)\mathop{}\!\mathrm{d}\tau$, which leads to the relation $\mathop{}\!\mathrm{d}B(t)=\xi(t)\mathop{}\!\mathrm{d}t$ between their differentials in the limit $\Delta t\to 0$. Analogously, the increment of the time-changed Brownian motion can be related to $\xi$ by the equation: $B(S(t+\Delta t))-B(S(t))=\int_{S(t)}^{S(t+\Delta t)}\xi(\tau)\mathop{}\!\mathrm{d}\tau$. Moreover, recalling that the paths of $S$ are continuous and monotonically increasing, we find that $\Delta S(t)=S(t+\Delta t)-S(t)\to 0$ in the limit $\Delta t\to 0$, i.e., in such limit we obtain: $\mathop{}\!\mathrm{d}B(S(t))=\xi(S(t))\mathop{}\!\mathrm{d}S(t)=\xi(S(t))\dot{S}(t)\mathop{}\!\mathrm{d}t$.

As discussed earlier in this section, $Y$ can be shown to be a semimartingale, as long as the parent process $X$ in Eq. (15)(left) is a semimartingale jacod714calcul . In our specific case, $X$ is a Brownian diffusive process, i.e., it satisfies this property. Moreover, thanks to the continuity of the stochastic paths of $S$, both the process $Y$ and its general functional $W$, defined as in Eq. (1), have continuous trajectories kobayashi2011stochastic . Thus, the Itô formula of $Y$ is given by Eq. (44), where its quadratic variation can be computed by employing Eqs. (52a, 52b) kobayashi2011stochastic and recalling that for $X$ normal diffusive we have: $[X,X]_{t}=\int_{0}^{t}\sigma^{2}(X(s))\mathop{}\!\mathrm{d}s$. Thus, we can write the following:

[TABLE]

We note that this same result was also derived in magdziarz2010path with a different approach. Finally, Eq. (55) leads to the following equation for the infinitesimal increment of the quadratic variation of the process $Y$:

[TABLE]

IV Derivation of the generalized Feynman-Kac formula

With these mathematical preliminaries in place, we can derive the generalized FK equation for quantities $\widehat{P}$ of the form Eq. (2), with $W$ given as the general functional Eq. (1). The underlying stochastic process $Y$ is assumed to be an anomalous process with general waiting times, described by the coupled Langevin Eqs. (15), with the noise $\eta$ specified by its characteristic functional Eq. (16). For pedagogical reasons, we first consider the case of a purely space-dependent force $F(y)$ in Eq. (15)(left) and a time independent functional $U(x,t)=U(x)$ in Eq. (1). In this specific case, a brief discussion of this derivation has been presented previously in cairoli2015anomalous . Here, we provide the full details of this calculation. We then discuss its extension to space- and time-dependent forces (Sec. IV.2) and time-dependent functionals (Sec. IV.4). The latter case results in a set of coupled integro-differential evolution equations for $\widehat{P}$.

IV.1 Space-dependent forces

We here consider the case of a purely space-dependent force $F(y)$ in Eq. (15)(left) and a time independent functional $U(x,t)=U(x)$ in Eq. (1). We start from the two dimensional joint process $Z(t)=(Y(t),W(t))$. As suggested in Sec. III.3, the process $Z$ is a semi-martingale with continuous paths, as also $Y$ and $W$. Thus, its Itô formula (for a general smooth function $f$) is obtained by adapting Eq. (45) kunita1997stochastic and is given explicitly by:

[TABLE]

In order to simplify this equation we need the following ingredients: (i) the time-discretised form of Eq. (54) that expresses the increments of $Y$ in terms of the time-change increments $\mathop{}\!\mathrm{d}S(t)$; (ii) the differential increment of the quadratic variation of $Y$ in Eq. (56); (iii) the quadratic variation $[W,W]_{t}$ and covariation $[Y,W]_{t}$, which are both null as $W$ is a finite variation process (Appendix C). Further recalling from Eq. (1) that $\mathop{}\!\mathrm{d}W\!(t)=U(Y(t))\mathop{}\!\mathrm{d}t$, we obtain:

[TABLE]

The equation for the double Fourier transform of the joint PDF $\widehat{P}(p,k,t)$ can be derived by evaluating Eq. (58) for $f(Z(t))=e^{i\,k\,Y(t)+i\,p\,W(t)}$. Specifically, we obtain:

[TABLE]

Finally, we need to take the ensemble average over the realizations of both $\xi$ and $\eta$, the latter determining the realizations of the process $S$. Within the Itô prescription, the last integral in the rhs of Eq. (59) cancels out. This is briefly proven in the following. Let us introduce a finite time-discretisation with mesh $\Delta t$ and let $N=t/\Delta t$. We denote: $Z_{i}=Z(t_{i})$, $S_{i}=S(t_{i})$ and $Y_{i}=Y(t_{i})$. The stochastic integral can be written as:

[TABLE]

Let us take the average over $\xi$ first. For each fixed realization of $S$ we can then write: $\left\langle f(Z_{i})\sigma(Y_{i})[B(S_{i+1})-B(S_{i})]\right\rangle=\left\langle f(Z_{i})\sigma(Y_{i})\right\rangle\left\langle B(S_{i+1})-B(S_{i})\right\rangle=0$ which is due to (i) the independence of the increments of $B$, that enables us to factorize the average because both $Z_{i}$, $Y_{i}$ only depends on its previous increments, and (ii) to the null first moment of $B$. Thus, the averaged Eq. (59) reduces to the following:

[TABLE]

where in the second integral the Fourier transform of the FP operator of Eq. (15)(left): $\mathcal{L}_{\rm FP}(y)=-\frac{\partial}{\partial y}F(y)+\frac{1}{2}\frac{\partial^{2}}{\partial y^{2}}\sigma^{2}(y)$ appears. Further recalling that the inverse Fourier transform of $f(Z(t))$ is equal to $e^{i\,p\,W(t)}\delta(y-Y(t))$ and by using the properties of the delta function, we can derive from Eq. (61) the following equation for $\widehat{P}(p,y,t)$ of Eq. (2):

[TABLE]

To close the equation, we need to relate the averaged stochastic integral in Eq. (62) to $\widehat{P}(p,y,t)$. To this aim, we first write $W$ as a subordinated process with the change of variables $\tau=S(r)$, i.e., $T(\tau)=r$, in Eq. (1):

[TABLE]

where the noise $\eta(s)$ explicitly appears from Eq. (15)(right). Thus, by employing the property $1=\int_{0}^{+\infty}\delta(s-S(t))\mathop{}\!\mathrm{d}s$ and then considering the same discretisation scheme and notation used to derive Eq. (60), we obtain:

[TABLE]

where (i) the continuity of the paths of $S$ implies that no jump terms appear in the stochastic integral (Appendix B) and (ii) we used Eq. (49) to relate the stochastic increments of $S$ to those of $T$. If we take the average over the realizations of the two noises $\eta$ and $\xi$ of Eq. (64) and then the time derivative of the resulting expression, we obtain:

[TABLE]

Furthermore, the rhs side of Eq. (65) can be related in Laplace space to the joint PDF $\widehat{P}(p,y,t)$. By using again the representation of $W$ of Eqs. (63) and the property used to derive Eq. (64), $\widehat{P}(p,y,t)$ can be rewritten as follows:

[TABLE]

Written in this form, its Laplace transform can be computed straightforwardly. Indeed, recalling Eqs. (15, 47), we can derive the Laplace transform of $\delta(s-S(t))$ as follows:

[TABLE]

such that the Laplace transform of Eq. (66) is given by:

[TABLE]

In Eq. (68), we explicitly highlighted that the ensemble average is made over the two different noises $\eta$ and $\xi$, whose independence allows us to change arbitrarily the order in which such averages are performed. This flexibility can be readily employed to simplify Eq. (68) by expressing the $\eta(s)$-dependent part of the integrand as a derivative of the characteristic functional $G[k(l)]$ in Eq. (16). Indeed, by performing the average with respect to $\eta(s)$ first and recalling that $X$ does not depend on it, this implying that the delta function can be taken out of such average, the only quantity needed to be computed is $\left\langle\eta(s)\,\exp{[-\lambda\,T(s)+i\,p\,A(s)]}\right\rangle$. This can be obtained as follows:

[TABLE]

where we used the characteristic functional Eq. (16) with the test function $k(l)=\left[\lambda-i\,p\,U(X(l))\right]\Theta(s-l)$. Substituting Eq. (69) back into Eq. (68), we derive the following relation:

[TABLE]

where now the brackets denote again an average over both $\eta$ and $\xi$. We note that the Laplace transform of the rhs of Eq. (65) is equal to the integral of Eq. (70). Thus, by expressing it in terms of ${\widehat{\widetilde{P}\mkern 3.5mu}\mkern-3.5mu}{}(p,y,\lambda)$, taking its inverse Laplace transform and substituting it back in Eq. (62), we derive the generalized FK formula:

[TABLE]

where the memory kernel is related to $\Phi$ trough the following relation (in Laplace space):

[TABLE]

Eq. (71) highlights that the non-Markovian features of the underlying anomalous process $Y$ result in a temporal memory that is directly related to both the statistics of the waiting times, expressed by the Laplace exponent $\Phi$, and the $y$-coordinate via the function $U$. Consequently, the integral operator expressing the temporal memory does not commute with the Fokker-Planck operator. In the specific case of $\Phi$ as in Eq. (41), i.e., $T$ is a one-sided tempered Lévy stable process, Eq. (71) has also recently been confirmed by using a master equation approach wu2016tempered .

IV.2 Space-time dependent forces

In the presence of both space- and time-dependent forces and no multiplicative term, Eq. (71) with the substitution $\mathcal{L}_{\rm FP}(y)\to\mathcal{L}_{\rm FP}(y,t)=-\frac{\partial}{\partial y}F(y,t)+\frac{\sigma}{2}\frac{\partial^{2}}{\partial y^{2}}$ was already proved in the specific case of CTRWs with power-law waiting times starting from a master equation approach in carmi2011fractional . As discussed in Sec. II, the time dependence in the external force is introduced by making $F$ depend explicitly on $T$, i.e., we consider the general dynamics described by the subordinated Langevin Eqs. (15). There are two main differences with the time independent case. On the one hand, the processes $X$ and $T$ are no longer independent, such that the previous derivation of Eq. (70) does not hold any more. Specifically, the delta function in Eq. (68) needs to be kept inside the average over the realizations of $\eta$. On the other hand, while in the time independent case the stochastic paths of $X$, and consequently those of $F$, have continuous paths, in the case of Eq. (15)(left), due to the explicit dependence of $X$ on the Lévy process $T$, both its paths and those of $F$ are generally càdlàg, with random jumps occurring in correspondence to those of $T$. Nevertheless, thanks to the finite variation of $T$, both $X$ and the time-changed process $Y$ are still semimartingales. Thus, we can integrate Eq. (15)(left) as below:

[TABLE]

We remark that the integral over $F$ is done with respect to a process with finite variation and continuous paths, i.e. a deterministic drift (Lebesgue measure), such that the contribution from the random jumps of $T$ is still null. As in the time independent case discussed earlier, we can use directly the time-change to write an equation for $Y$:

[TABLE]

where we employed again Eq. (52a). After taking its time derivative, we derive the following equation kobayashi2011stochastic :

[TABLE]

This result elucidates that $Y$ has still continuous paths and that both Eqs. (44, 56) still hold. Similar arguments as in the previous derivation for the time independent case can be made, leading to the same Eq. (62) with $\mathcal{L}_{\rm FP}(y,t)$ and the same averaged stochastic integral, which needs to be related to the joint PDF. As already highlighted, the proof of Eq. (70) needs a more detailed analysis, as both $X$ and $W$ now depend on the realizations of $\eta$. Starting from Eq. (68), we can first rewrite the explicit $\eta(s)$ dependence as a time derivative of the exponential function as follows:

[TABLE]

where we use again Eq. (63) and we employ the properties of the delta function to factorize the term $[\lambda-i\,p\,U(y)]^{-1}$ out of the integral. Differently from the time independent case, the factors inside the ensemble average can no longer be separated, i.e., we cannot compute directly such term by means of Eq. (16). Nevertheless, such expression can still be simplified if we look at its discretised form. We consider a partition $\pi=\{0=s_{0}<s_{1}<\ldots<s_{n}=s\}$ of the interval $[0,s]$ with constant mesh $\Delta s$ and $n=s/\Delta s$. We denote: $X(s_{i})=X_{i}$ and $\eta(s_{i})=\eta_{i}$. We recall that $\eta_{i}$ are RVs with characteristic function specified by $\Phi$, such that $\Delta T_{i}=\eta_{i}\,\Delta s$ is the corresponding $T$ increment (according to Eq. (15)(right)). As the delta function in Eq. (76) only imposes a condition on the final point, we can write:

[TABLE]

where in the first line we discretise the derivative in the operational time $s$ and in the third one we factorized the average over the last increment $\Delta T=\Delta s\,\eta_{n+1}$. This is allowed because (i) $U(X_{n-1})$ only depends on the increments of the process $T$ up to $n$, which are independent on the RV $\eta_{n+1}$, and (ii) the end-point value is conditioned to $y$, i.e. $U(X_{n})=U(y)$, which is no longer a RV. The average is then computed with Eq. (16):

[TABLE]

Substituting this term into Eq. (77) and taking the continuum limit $\Delta s\to 0$, we obtain:

[TABLE]

leading with Eq. (76) to the same relation Eq. (70) also in the case of both space- and time-dependent forces. The rest of the derivation follows as in the time independent case. Thus, we have shown that Eq. (71) with the substitution $\mathcal{L}_{\rm FP}(y)\to\mathcal{L}_{\rm FP}(y,t)=-\frac{\partial}{\partial y}F(y,t)+\frac{1}{2}\frac{\partial^{2}}{\partial y^{2}}\sigma^{2}(y)$ is the generalized FK formula of processes described by the subordinated Langevin Eqs. (15), where external forces are allowed to depend on time, as well as on space. Thus, our formalism naturally provides a solution to the issue of the position of the FP operator with respect to the memory integral in a more general framework than CTRWs, in which case it has long been debated heinsalu2007use ; magdziarz2008equivalence ; weron2008modeling ; henry2010fractional .

IV.3 Special Cases and Extensions

Our proposed FK Equation (71) recovers several different equations earlier derived in the literature for specific choices of the waiting time distribution and/or of the function $U$. We summarize these special cases below.

The generalized Fokker-Planck Equation. If we set $p=0$, we find a generalized Fokker-Planck equation for the position PDF $f_{Y}$ magdziarz2009langevin :

[TABLE] 2. 2.

The generalized Klein-Kramers Equation. If we set $U(x)=x$ in Eq. (1), $Y$ and $W$ correspond respectively to the velocity and the position of an anomalous diffusing particle. Thus, after inverse Fourier transform, Eq. (71) yields a generalized fractional Klein-Kramers equation, which extends the result of Friedrich2006Anomalous ; Friedrich2006Exact :

[TABLE]

The shift of the position sample variable in the memory integral elucidates the presence of the same retardation effects of Friedrich2006Anomalous ; Friedrich2006Exact . In particular, our derivation highlights that the stochastic dynamics underlying Eq. (81) is given by the coupled Langevin Eqs. (15), which has been conjectured without proof in eule2007langevin . 3. 3.

Normal Diffusion. This case is obtained with $\Phi(\lambda)=\lambda$, i.e., $K(t)=1$. In this case, Eqs. (71, 80, 81) reduce respectively to the ordinary FK , Fokker-Planck and Klein-Kramers equation risken1989fokker . 4. 4.

CTRWs with power-law waiting times. This case is obtained by setting $\Phi(\lambda)=\lambda^{\alpha}$ with $0<\alpha<1$ [Eq. (39)]. In this case, $K(t)=t^{\alpha-1}/\Gamma(\alpha)$, such that the integral operator in Eq. (71) specifies to

[TABLE]

which is the fractional substantial derivative introduced in Friedrich2006Anomalous ; Friedrich2006Exact ; turgeman2009fractional ; carmi2011fractional ; Orzel2011FKKE . With such choice, Eqs. (71, 81) become the fractional FK equation turgeman2009fractional ; carmi2011fractional ; Orzel2011FKKE and the fractional Klein-Kramers equation Friedrich2006Anomalous ; Friedrich2006Exact respectively. If we set $p=0$, the previous operator further reduces to the Riemann-Liouville fractional derivative, i.e., Eq. (80) becomes the fractional diffusion equation metzler1999deriving ; magdziarz2008equivalence . 5. 5.

CTRWs with tempered Lévy-stable distributed waiting times. This case is obtained by setting $\Phi(\lambda)=(\lambda+\mu)^{\alpha}-\mu^{\alpha}$ with $\alpha$ as in the previous case and the tempering index $\mu\in\mathbb{R}^{+}$ [Eq. (41)]. The corresponding memory kernel in Eqs. (71, 80) is $K(t)=e^{-\mu t}t^{\alpha-1}E_{\alpha,\alpha}[(\mu t)^{\alpha}]$, with $E_{\alpha,\alpha}$ being a two-parameter Mittag-Leffler function. This specific case has also been recently discussed in wu2016tempered by solving directly for the Laplace-Fourier transform of the joint PDF ${\widehat{\widetilde{P}\mkern 3.5mu}\mkern-3.5mu}{}(p,k,\lambda)$ of a suitable CTRW and then taking its diffusive limit, i.e., $(k,\lambda)\to(0,0)$. Here, we prove the equivalence of our own result and the approach therein by deriving such a limit solution. For simplicity, we restrict to time-independent external forces. Let us take the Fourier-Laplace transform of Eq. (71). Recalling that the functions $U$ and $F$ are smooth and using the convolution theorem $yf(y)\to-i\frac{\partial}{\partial k}\widehat{f}(k)$, their Fourier transforms can be expressed as $U(y)\to U\!\left(-i\frac{\partial}{\partial k}\right)$ and $F(y)\to F\!\left(-i\frac{\partial}{\partial k}\right)$. This is understood by first Taylor expanding these functions, then Fourier transforming each term separately and finally by re-summing the series expansions. Thus, assuming the initial condition $P(w,y,t=0)=\delta(w)\delta(y)$, we obtain from Eq. (71)

[TABLE]

Rearranging the terms, we can rewrite it as

[TABLE]

We note that the terms in front of ${\widehat{\widetilde{P}\mkern 3.5mu}\mkern-3.5mu}{}$ are operators in the Fourier variable $k$, that do not commute in general. Therefore, applying their inverse to both sides of the previous equation in the correct order, we derive:

[TABLE]

Substituting the $\Phi$ prescribed, we recover Eq. (3) (for $F=0$) and Eq. (16) of wu2016tempered . Eq. (85) is the formal solution in the diffusive limit of the joint PDF of a CTRW with waiting time distribution $\widetilde{\psi}(\lambda)=e^{-\Phi(\lambda)}$ turgeman2009fractional ; carmi2011fractional . 6. 6.

Multiplicative $X-$process with general $\alpha$-prescription. We consider the set of subordinated Langevin equations:

[TABLE]

where $\star$ denotes a generalized $\alpha$ prescription in the definition of the stochastic integral (as in Eq. (112a) in Appendix A). However, the resulting process $X$ is equivalently described by Eqs. (15), i.e., with the ordinary Itô prescription, by using the mapping given by Eqs.(120a, 120b) (details are presented in Appendix A). Thus, Eqs. (71, 72) still hold for the subordinated Eqs. (86) with the modified FP operator:

[TABLE]

The dynamics of $W$ in Eq. (1), when $Y$ is obtained by subordination of a process $X$ of the type described by Eq. (86), exhibits peculiar behavior, e.g. Lévy flight dynamics lubashevsky2009realization ; lubashevsky2009continuous , already in the Brownian limit, i.e., $\Phi(\lambda)=\lambda$. This motivates our interest in extending our generalized FK Eq. (71) to such types of $X$-processes.

IV.4 Explicit time dependence in the functional

As highlighted previously in Sec. I a generic time-dependent protocol driving a system out of equilibrium leads both to a time-dependent force and to an explicit time-dependence in the functional, which defines the accumulated mechanical work. To the extent of our knowledge, such time-dependent functionals in the form of Eq. (1) have so far not been discussed in the literature of anomalous diffusive processes. We here address this issue by considering $W(t)$ given by Eq. (1), where the dynamics of the underlying process $Y$ is represented by the subordinated Langevin Eqs. (15). The time-dependence in $U$ does not modify the properties of $W$, i.e., it still has finite variation and continuous paths. Thus, the Itô formula for the two-dimensional semimartingale $Z(t)=(Y(t),W(t))$ is given by Eq. (58) with the substitution $U(Y(\tau))\rightarrow U(Y(\tau),\tau)$. Consequently, by repeating a similar calculation as that presented in Sec. IV, we obtain the following equation:

[TABLE]

By using the same change of variables employed for Eq. (63), one can rewrite $W$ as below:

[TABLE]

To proceed, we note that (i) Eq. (64) does not involve the explicit definition of $A$, such that the fundamental Eq. (65) also holds in this case, and (ii) the double average in Eq. (68) can no longer be factorized, because of the dependence of $A$ on the process $T$. Nevertheless, we can address this issue with an argument similar to that presented in Sec. IV.2. Thus, we first modify Eq. (68) as follows:

[TABLE]

Secondly, we manipulate the ensemble average appearing in its rhs by considering its discretised form. By employing the same notation used to derive Eq. (77), we obtain:

[TABLE]

where we explicitly separate the average over the last increment $\eta_{N+1}$ and that over the increments of both $\xi_{j}$ and $\eta_{j}$ for $j\leq n$. Differently from Eq. (77), these two averages cannot be factorized, because the RV $T_{n}$ depends on all the increments $\eta_{j}$. However, as $\eta_{n+1}$ is independent on $\eta_{j}$, the internal average can be solved by using Eq. (16):

[TABLE]

such that, by substituting it into Eq. (91) and taking the continuum limit $\Delta s\to 0$, leads to the following relation:

[TABLE]

Finally, if we substitute it back into Eq. (90), we obtain:

[TABLE]

In order to close the evolution equation, we introduce the auxiliary function $\Omega$, which is defined in Laplace space as:

[TABLE]

such that the inverse Laplace transform (denoted as $\mathcal{L}^{-1}$) of the time dependent terms of Eq. (94) can be written as

[TABLE]

By employing this result in Eq. (94), we obtain the equation:

[TABLE]

where in the second line we changed the order of integration. Remarkably, the term in square brackets is the same integral in the rhs of Eq. (65). By using Eqs. (65, 97), we obtain that Eq. (88) is equivalent to the coupled equations:

[TABLE]

Here, $H$ is an auxiliary function that is coupled to the joint PDF $\widehat{P}$. We note that in the time independent case $U(x,t)=U(x)$, Eq. (99) reduces to a Laplace convolution. As a consequence, its Laplace transform factorizes, such that it can be solved explicitly for $H$, which can then be substituted in Eq. (98) to obtain a single closed equation, which recovers the generalized FK Eq. (71) with a space- and time-dependent force.

V Application to a non-equilibrium particle model

Since our framework includes a space- and time-dependent force as well as an explicit time-dependent functional, we can apply it to calculate the work fluctuations of an anomalous system driven by an arbitrary non-equilibrium protocol $q(t)$. Assuming a time-dependent potential $V(x,q(t))$, we consider the dynamics of Eqs. (15) with $F(x,t)=-\partial_{x}V(x,q(t))$ and $W(t)=\int_{0}^{t}\mathrm{d}r\,\partial_{q}V(Y(r),q(r))\dot{q}(r)$. One of the simplest examples is that of a potential moving at constant velocity $v_{0}$. In this case we can set $q(t)=v_{0}t$ and $V(x,q(t))=V(x-q(t))$. This system represents the simplest pure out-of-equilibrium model, where a steady-state can be reached from the balance between dissipative forces, i.e., the friction of the surrounding fluid, and driving forces, i.e., the time dependent force due to the moving potential. Thus, it provides an easily solvable setup, where the applicability of fluctuation theorems for the accumulated mechanical work done by the system, in particular the so-called steady-state fluctuation theorem evans1993probability ; gallavotti1995dynamical ; gallavotti1995dynamical ; gallavotti1995reversible ; kurchan1998fluctuation ; lebowitz1999gallavotti ; seifert2012stochastic , can be tested. A widely studied case is that of an harmonic potential $V(x,t)=\gamma\,(x-v_{0}\,t)^{2}/2$, such that we have in Eqs. (15)

[TABLE]

We also set $\sigma(x)=\sqrt{\sigma}$, where $\sigma$ is a positive real constant. In the limit case of normal diffusive dynamics, this model has been extensively studied both theoretically and experimentally wang2002experimental ; van2003stationary ; trepagnier2004experimental ; taniguchi2007onsager ; taniguchi2008nonequilibrium ; gomez2010steady ; aquino2013power ; mestres2014realization . The work fluctuations in the steady-state regime are described by the large deviation function $I(w)=-\lim_{t\to\infty}t^{-1}P(w,t)$, which has been calculated explicitly not only in the normal diffusive regime, but also when the random force exerted by the external bath is described by either a Lévy or a Poisson shot noise touchette2007fluctuation ; touchette2009anomalous ; baule2009steady . The solution of this paradigmatic model for anomalous dynamics of the type described in Sec. II has so far not been obtained.

In the anomalous case, the Eqs. (98, 99) can be applied in principle to calculate the joint PDF $\widehat{P}(p,y,t)$ [Eq. (2)], and consequently the large deviation function $I(w)$. Due to the linear form of $W(t)$ in Eq. (100)(right), we can also use the simpler FK Eq. (71), since the time-dependence in the functional can be separated as below:

[TABLE]

However, even for the linear dynamics of the dragged harmonic potential, arguably one of the simplest ways to impose a space- and time-dependent non-equilibrium drive, the linear functional case could not be fully solved so far.

Here, we study the first and second moment of both the position and the work. We compute the first two moments of $Y$, which will be needed later to compute the corresponding ones of $W$. The distribution of $Y$ is given by the generalized Fokker-Planck Equation (80), once we account for the correct time-dependent Fokker-Planck operator:

[TABLE]

with the memory kernel $K$ specified by Eq. (72). The general $n$-th order moment of $Y$ can be computed with the following standard procedure: (i) we take the Laplace transform of Eq. (102) with the remark that the linear time dependent term produces a derivative in the Laplace variable; (ii) we multiply both its sides by $y^{n}$ (clearly $n=1,2$ for the first and second moment respectively); (iii) we perform the ensemble average, i.e., we integrate in $y$ both sides of the resulting equation. We remark that in Laplace space one does not need to specify beforehand the waiting time distribution, i.e., we can derive results for general $\Phi$. For the first moment, we obtain the following formula:

[TABLE]

In the case of CTRWs, i.e., $\Phi(\lambda)=\lambda^{\alpha}$, Eq. (103) reduces to $\langle\widetilde{Y}(\lambda)\rangle=\left[y_{0}+\gamma\,v_{0}\,\alpha\,\lambda^{-1-\alpha}\right]/[\lambda+\gamma\,\lambda^{1-\alpha}]$, whose inverse Laplace transform can be computed analytically: $\langle Y(t)\rangle=y_{0}\,E_{\alpha}(-\gamma\,t^{\alpha})+v_{0}\,\alpha\,t\,[1-E_{\alpha,2}\!\left(-\gamma\,t^{\alpha}\right)]$. Furthermore, if we set $\alpha=1$ (Brownian limit) and recall that in such limit $t\,E_{\alpha,2}\!\left(-\gamma\,t^{\alpha}\right)\rightarrow(1-e^{-\gamma\,t})/\gamma$, we obtain the expected result: $\langle Y(t)\rangle=y_{0}\,e^{-\gamma\,t}+v_{0}\,t-[v_{0}(1-e^{-\gamma\,t})]/\gamma$ risken1989fokker ; gardiner1985handbook . For the second order moment, we need to compute the quantity: $\partial_{\lambda}[\lambda\,\langle\widetilde{Y}(\lambda)\rangle/\Phi(\lambda)]$, which is due to the time dependent force term. By using Eq. (103), we find:

[TABLE]

In the Lévy-stable case, $\Phi(\lambda)=\lambda^{\alpha}$, Eqs. (104a, 104b) can be shown to reduce to the following ones:

[TABLE]

which can be Laplace inverse transformed analytically as below:

[TABLE]

As a sanity check, by setting $\alpha=1$ in Eq. (106) we obtain: $\langle Y^{2}(t)\rangle=y_{0}^{2}\,e^{-2\,\gamma\,t}+v_{0}^{2}\,[(1-e^{-\gamma\,t})^{2}/\gamma^{2}]+v_{0}^{2}\,t^{2}+2\,v_{0}\,y_{0}\,t\,e^{-\gamma\,t}-2\,v_{0}\,y_{0}\,e^{-\gamma\,t}\,[(1-e^{-\gamma\,t})/\gamma]-2\,v_{0}^{2}\,t\,[(1-e^{-\gamma\,t})/\gamma]+2\,\sigma\,[(1-e^{-2\,\gamma\,t})/\gamma]$, which is the expected Brownian limit. In the case of the mechanical work, the first two moments are obtained by using Eq. (101). Thus, we have:

[TABLE]

where both $\langle Q(t)\rangle$ and $\langle Q^{2}(t)\rangle$ can be computed analytically. On the one hand, thanks to the linearity of the functional, we find: $\langle\widetilde{Q}(\lambda)\rangle=\langle\widetilde{Y}(\lambda)\rangle/\lambda$, such that Eq. (103) can be employed to derive a closed analytic expression. On the other hand, the second order moment of $Q$ is derived by exploiting the following FK equation (Eq. (71) adapted explicitly to the case considered here):

[TABLE]

To this aim, we need a procedure to compute joint moments of the type $\langle[Y(t)]^{m}\,[Q(t)]^{n}\rangle$. Recalling that the Fourier transform of the joint PDF of $Q$ and $Y$ is equal to $\widehat{P}(p,y,t)=\langle e^{i\,p\,Q(t)}\,\delta(y-Y(t))\rangle$, the general $n$-th order moment of $Q$ is given by $\langle[Q(t)]^{n}\rangle=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}a^{n}P(a,y,t)\mathop{}\!\mathrm{d}a\mathop{}\!\mathrm{d}y=(-i\,\partial_{p})^{n}\int_{-\infty}^{+\infty}\widehat{P}(p,y,t)\mathop{}\!\mathrm{d}y|_{p=0}$. In the case of the joint moment, one simply needs to include a factor $y^{m}$ in the integral, i.e., $\langle[Y(t)]^{m}[Q(t)]^{n}\rangle=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}a^{n}\,y^{m}P(a,y,t)\mathop{}\!\mathrm{d}a\mathop{}\!\mathrm{d}y=(-i\,\partial_{p})^{n}\int_{-\infty}^{+\infty}y^{m}\,\widehat{P}(p,y,t)\mathop{}\!\mathrm{d}y|_{p=0}$. Thus, the strategy to compute such moments is the following: (i) we take the Laplace transform of Eq. (108) in order to express the integral term by means of $\Phi$; (ii) we multiply each side of the resulting equation by the corresponding power of $y$ and take its integral; (iii) we make the corresponding derivative in $p$ and evaluate the expression for $p=0$. In the specific case of the second order moment of $Q$, we obtain:

[TABLE]

where we define the following auxiliary functions:

[TABLE]

Our analytic results for the $Y$ and $W$ moments are in full agreement with numerical simulations, as shown in Fig. 2, where we choose the waiting time process to be tempered Lévy-stable [$\Phi$ as in Eq. (41)]. The normal diffusive regime corresponds to a plateau in both Fig. 2a and 2b, since both moments are plotted rescaled by $1/t^{2}$ and become ballistic for large times van2003stationary . Interestingly, we observe a crossover scaling for all values of $\mu$ including the pure Lévy-stable regime ($\mu=0$). For short times, one generally observes the subdiffusive scaling $\langle Y^{2}(t)\rangle\sim t^{\alpha}$ for the position coordinate and the superdiffusive scaling $\langle W^{2}(t)\rangle\sim t^{\alpha+2}$ for the work. For long times $\langle Y^{2}(t)\rangle$ converges to the ballistic scaling of the normal diffusive scenario for all $\mu$, though with different ballistic diffusion coefficient between the CTRW case $\mu=0$ and that of finite $\mu$. The behavior of $W$ is instead qualitatively different, as we find $\langle W^{2}(t)\rangle\sim t^{4}$ for $\mu=0$ and $\langle W^{2}(t)\rangle\sim t^{2}$ for finite $\mu$. Different intermediate scaling are observed depending on $\mu$.

VI Conclusions and open questions

Despite the widespread occurrence of anomalous diffusive processes in physical, chemical and biological systems, only recently the properties of their general observables have been investigated through the systematic derivation of a FK-type equation. In this paper we demonstrated how such an equation can be derived from the stochastic description of a CTRW with generalized waiting times in the diffusive limit. Thus, our results extend the correspondence between Langevin equations and deterministic partial differential equations of the original Feynman-Kac theorem to the anomalous regime. This has been obtained by expressing the waiting time statistics in terms of a general one-sided Lévy process, whose Laplace exponent turns out to be directly related to the memory kernel of the generalized FK equation. Formally, the CTRW with generalized waiting times is expressed as a normal diffusive process, subordinated by the inverse of the Lévy process, which is itself generally not a Lévy process, but a semimartingale. Consequently, our derivation requires a basic knowledge of Lévy processes, semimartingales and their stochastic calculus, which we here provided in a pedagogical way.

While the case of space- and time-dependent forces is a non-trivial extension of our previous results in Ref. cairoli2015anomalous , it is manifest only in a modified Fokker-Planck operator in the generalized FK equation as expected. On the other hand, the case of an explicit time dependence in the functional leads to a new class of coupled integro-differential evolution equations for the joint PDF in Fourier space. The main challenge now is to develop methods to derive analytical solutions of such equations, that could be applied to physically relevant situations. In fact, even for the simplest scenario of time independent forces and functionals, explicit solutions are sparse and have been restricted to moments or asymptotic expressions for the simplest observables with underlying Lévy-stable Friedrich2006Anomalous ; Friedrich2006Exact ; turgeman2009fractional ; carmi2011fractional ; carmi2010distributions and tempered Lévy-stable waiting time processes cairoli2015anomalous ; wu2016tempered .

Open questions on a conceptual level concern the derivation of backward FK-type equations in our framework and the inclusion of a time dependence in the multiplicative diffusion term. In the master equation approach, backward FK equations, i.e., where the spatial derivatives act on the space coordinate at the initial time, are straightforward to derive and are particularly relevant for occupation time problems turgeman2009fractional ; carmi2011fractional ; carmi2010distributions ; wu2016tempered . On the contrary, in the subordination framework, the backward equations are much more challenging to treat. Conversely, the case of a time- dependent diffusion term can be treated along the lines presented here, but details are left for future work.

Our results are in particular applicable to the stochastic thermodynamics of anomalous processes. Previous studies of fluctuation theorems in anomalous subdiffusive systems focused on work induced by a constant force such that the work statistics are equivalent to that of the spatial coordinate itself chechkin2009fluctuation ; chechkin2012normal ; dieterich2015fluctuation . On the other hand, the mechanical work imposed by a non-equilibrium driving in a generic situation is naturally captured by our framework. A detailed discussion of work fluctuations and the associated fluctuation theorems relies on a knowledge of the large deviation function, which is currently out of reach already for the simple model of an anomalous particle in a moving harmonic potential here discussed. A further study of this paradigmatic model is certainly valuable to gain fundamental insight into the interplay of non-equilibrium driving and complex waiting time processes. Moreover, such a system can be implemented in a straightforward way in experiments, e.g., by immersing a tracer particle in a complex fluid environment and dragging it with optical tweezers tassieri2016microrheology . Thus, our results pave the way for the theoretical investigation of non-equilibrium processes of this type.

APPENDIX A Relation between generalized $\bm{\alpha}$ and Itô Prescription

We review the definition of the stochastic integral with respect to a Brownian motion $B(t)$. Let us introduce (i) a process $Y$ a.s. continuous, (ii) a Brownian motion $B$ on the time interval $[0,t]$ and (iii) a partition $\pi=\{0=t_{0}<t_{1}<\ldots<t_{N}=t\}$ of the interval $[0,t]$ with finite mesh $|\pi|$, such that $|\pi|\to 0$. Thus, we can define the stochastic integral of $Y$ with respect to the increments of $B$ the following stochastic process revuz1999continuous ; applebaum2009levy ; karatzas2012brownian :

[TABLE]

However, the choice of the specific time at which we evaluate the integrand process $Y$ in Eq. (111) is arbitrarily chosen. There, this is the earlier time $t_{i}$. This specific choice is called the Itô prescription, but in general one can choose any point in the interval $[t_{i},t_{i+1}]$. Each of these different choices generate integrals with completely different properties. A general definition of the stochastic integral accounting for all the different prescriptions is given in terms of a parameter $\alpha\in[0,1]$ as follows lau2007state ; lubashevsky2009realization :

[TABLE]

For $\alpha=0$, we recover the Itô prescription, whereas for $\alpha=1/2$ we obtain the Stratonovich prescription. The case $\alpha=1$ has also been discussed in hanggi1982stochastic ; klimontovich1990ito . Processes of this type will be denoted with $\star$.

We now show that 1D stochastic processes with general $\alpha$-prescription can be mapped into Itô processes by suitably choosing the coefficients of the Langevin equation. Specifically, we consider a process $Y(t)$ described by

[TABLE]

where we use respectively the generalized $\alpha$-prescription ($\alpha\neq 0$) as in Eq. (112b) or the Itô one. Our aim is to find suitable functions $a$, $b$, such that the integrated process is the same. We consider the integrated version of Eq. (113):

[TABLE]

where the stochastic integral is defined as in Eqs. (112a, 112b). Here, we used the relation $\mathop{}\!\mathrm{d}B(t)=\xi(t)\mathop{}\!\mathrm{d}t$ between the increments of a Brownian motion and the white Gaussian noise $\xi$ gardiner1985handbook . Our first task is to represent this term as an Itô stochastic integral. To this aim, let us consider the same partition $\pi$ as before and rewrite the auxiliary variable $Z$ as $Z(t_{i})=Y(t_{i})+\alpha\,[Y(t_{i+1})-Y(t_{i})]$. Thus, we can write:

[TABLE]

We note that $Z(t_{i})$ depends on $\Delta Y(t_{i})=Y(t_{i+1})-Y(t_{i})$, which can be expressed as an Itô increment by using the discretised version of Eq. (114), i.e., we find $\Delta Y(t_{i})=a(t_{i},Y(t_{i}))\Delta t_{i}+b(t_{i},Y(t_{i}))\Delta B(t_{i})$. To simplify the notation, we denote: $\Delta B(t_{i})=[B(t_{i+1})-B(t_{i})]$ and $\Delta t_{i}=(t_{i+1}-t_{i})$. Thus, we can employ such relation to express $\sigma$ as

[TABLE]

This result needs to be substituted back into Eq. (116). We can then further simplify such expression by recalling that $\left\langle\Delta B(t_{i})^{2}\right\rangle=\Delta t_{i}$, as they are Gaussian distributed by definition, and that the $\Delta t_{i}$ dependent term cancels out in the limit of null mesh. Thus we obtain:

[TABLE]

By using Eqs. (115, 118) together, we obtain:

[TABLE]

It is now clear that the mapping between the two processes is realised if we set:

[TABLE]

APPENDIX B Finite variation processes

We review definition and properties of continuous stochastic processes with paths of finite variation. Specifically, we will provide the definition of their stochastic integral and their Itô formula. As the time-change process $S$ defined in Sec. II belongs to this class, such notions are employed in the derivation of the generalized fractional FK Eq. (71).

As a preliminary step, we define the total variation of a real-valued function $g$ with support on an interval $[s,t]$. Thus, we consider a partition of the interval $\pi=\{s=t_{0}<t_{1}<\ldots<t_{N}=t\}$, whose mesh is given by the maximum of the lengths of the subintervals: $|\pi|=\max_{i=1,\ldots,N}|t_{i}-t_{i-1}|$ and compute the quantity:

[TABLE]

whose value clearly depends on the specific $\pi$ chosen. Let us now consider the set of all possible partitions $\mathcal{P}=\{\pi_{n}\}$ and the corresponding variations of $g$ with respect to them $\{V_{t}^{\pi_{n}}(g)\}$. The total variation of $g$ on $[s,t]$ is obtained by taking the supremum of this set:

[TABLE]

Thus, if $V_{t}(g)<\infty$, then $g$ is said to be of finite variation and $V_{t}(g)$ is the total variation of $g$ on the chosen interval; otherwise, it is said to have infinite variation. If $g$ is defined over all $\mathbb{R}$, then $g$ has finite variation if it is of finite variation on all closed intervals of $\mathbb{R}$. Clearly, if $g$ is a non decreasing function, then it is of finite variation, as $V_{t}(g)=g(t_{N})-g(t_{0})$. Conversely, if $g$ is of finite variation, we can always find two auxiliary non decreasing functions $g_{1}$ and $g_{2}$, such that $g=g_{1}+g_{2}$.

In a similar way, a stochastic process $Y$ is said to be of finite variation if its stochastic trajectories $Y(t)$ have finite variation almost surely, i.e., for each of its different realizations. An analogous definition holds in the opposite case of a process of infinite variation. We note that ordinary integrals (in Lebesgue sense) of a continuous stochastic process are also of finite variation.

Stochastic integrals with respect to these processes can be defined straightforwardly as Lebesgue-Stieltjes integral with the proper measure associated to $Y$, which exists due to the finite variation of their paths ash2000probability . In terms of Riemann sums, considering the same partition $\pi$ as before, the stochastic integral of an arbitrary function $H(t)$ with respect to $Y$ is defined as

[TABLE]

We note that $H$ does not need to be continuous, but its paths are required to be right continuous with left limits (càdlàg). We further remark that the stochastic integral in Eq. (123) can also be defined when $Y$ has general càdlàg paths, i.e., not continuous. However, in such case jump terms need to be properly accounted for. As this is not the case of $S$, we will not discuss it in this context. For a general differentiable function $f$ of $Y$ the Itô formula is

[TABLE]

This follows straightforwardly by considering again the partition $\pi$ and employing the mean value theorem:

[TABLE]

APPENDIX C The quadratic variation

We here define the quadratic variation of a process $Y(t)$ for $t\geq 0$ on a time interval $[s,t]$. Let us consider a partition of such interval $\pi=\{s=t_{0}<t_{1}<\ldots<t_{N}=t\}$ of mesh $|\pi|=\max_{i=1,\ldots,N}|t_{i}-t_{i-1}|$. Associated to $\pi$, we can define the following process:

[TABLE]

Clearly, $[Y,Y]^{\pi}_{t}$ depends both on the specific realization of $Y$ and on the partition chosen. To avoid this latter dependence, we study its properties in the limit $|\pi|\to 0$. Let us now consider sequences of partitions $\{\pi_{n}\}$, such that $|\pi_{n}|\to 0$, and compute the corresponding sequences $\{[Y,Y]^{\pi_{n}}_{t}\}$. If for every $t$ this latter sequence converges in probability to a finite value $[Y,Y]_{t}$ independent on the specific choice of $\{\pi_{n}\}$ a.s., then $[Y,Y]_{t}$ is a well-defined process called the quadratic variation of $Y$. As an example, we show that the quadratic variation of a process $Y$ with paths of finite variation exists and it is null. From Eq. (126) and for a given $\pi$, we can write:

[TABLE]

where $V_{t}^{\pi}(Y)$ is the variation of $Y$ as defined in Eq. (126), which is finite by assumption. Thus, the rhs converges to zero a.s. in the limit $|\pi|\to 0$. As this result holds independently of the specific $\pi$, we conclude that a.s. $[Y,Y]_{t}=0$.