Feynman-Kac equation for anomalous processes with space- and time-dependent forces Andrea Cairoli∗ and Adrian Baule† School of Mathematical Sciences, Queen Mary, University of London, Mile End Road, E1 4NS, UK (Dated: January 9, 2017) Functionals of a stochastic process Y(t) model many physical time-extensive observables, for in- stance particle positions, local and occupation times or accumulated mechanical work. When Y(t) isanormaldiffusiveprocess,theirstatisticsareobtainedasthesolutionofthecelebratedFeynman- 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) 7 1 isnon-Brownian,e.g.,ananomalousdiffusiveprocess,generalizationsoftheFeynman-Kacequation 0 that incorporate power-law or more general waiting time distributions of the underlying random 2 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 directly related to the memory kernel appear- n ing in the generalized Feynman-Kac equation. The corresponding anomalous processes have been a showntocapturenonlinearmeansquaredisplacementsexhibitingcrossoversbetweendifferentscal- J ingregimes,whichhavebeenobservedinnumerousexperimentsonbiologicalsystemslikemigrating 6 cellsordiffusingmacromoleculesinintracellularenvironments. However,thecasewherebothspace- and time-dependent forces drive the dynamics of the generalized anomalous process has not been ] h solved yet. Here, we present the missing derivation of the Feynman-Kac equation in such general c case by using the subordination technique. Furthermore, we discuss its extension to functionals e explicitly depending on time, which are of particular relevance for the stochastic thermodynamics m of anomalous diffusive systems. Exact results on the work fluctuations of a simple non-equilibrium modelareobtained. AnadditionalaimofthispaperistoprovideapedagogicalintroductiontoL´evy - t processes,semimartingalesandtheirassociatedstochasticcalculus,whichunderliethemathematical a formulation of anomalous diffusion as a subordinated process. t s . t a m - I. INTRODUCTION d n In experimental applications one typically measures physical observables W, whose time evolution is determined o c 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: 1 (cid:90) t v W(t)= U(Y(r),r)dr, (1) 1 0 4 where U(x,t) is some prescribed arbitrary function. If Y is a normal diffusive process, these functionals have been 6 employed to model many different physical phenomena by choosing the function U suitably either with an explicit 1 0 timedependenceorwithoutit. Forinstance, inthelinearcaseU(x,t)=x, withY interpretedasaparticle’svelocity, . W represents its position and the Y–W representation simply describes the stochastic evolution of the system in the 1 phase space [1]. If instead we choose U(x,t) = δ(x) and U(x,t) = Θ(x), W stands for the local and occupation 0 7 timerespectively[2–13]. OtherrelevantchoicesareU(x,t)=x2 withEq.(1)interpretedasaone-dimensionalspatial 1 integral,inwhichcaseW isinterpretedasthevarianceofafluctuatinginterface[14],andU(x,t)=e−βx withβ areal : positive parameter, which describes the dynamics of integrated stock prices of the Black-Scholes type [15]. Another v important class of functionals has been introduced in the context of the stochastic thermodynamics of driven small i X scale systems. Specifically, one is interested in the statistical properties of the accumulated mechanical work done r by the system when a non-equilibrium driving is imposed by a time-dependence in a potential V(x,q(t)) via some a prescribedtime-dependentprotocolq(t). Insuchascenario,assumingY tobethepositioncoordinate,themechanical work is defined by Eq. (1) with U(x,t)=∂ V(x,q(t))q˙(t) [16–20]. q In order to obtain the probability density function (PDF) of W one usually considers quantities of the form: (cid:68) (cid:69) P(cid:98)(p,y,t)= eipW(t)δ(y−Y(t)) (2) ∗ Presentaddress: DepartmentofBioengineering,ImperialCollegeLondon,SouthKensingtonCampus,SW72AZ,UK † Correspondingauthor: [email protected] 2 for a given initial condition Y(0)=y0. We note that P(cid:98) 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), P(cid:98) can be obtained by solving the Feynman-Kac (FK) equation [14]: ∂ ∂tP(cid:98)(p,y,t)=ipU(y,t)P(cid:98)(p,y,t)+LFP(y,t)P(cid:98)(p,y,t), (3) whereweintroducethemostgeneralFokker-Planckoperatorforaspace-timedependentforceanddiffusioncoefficient: L (y,t) = − ∂ F(y,t)+ 1 ∂2 σ2(y,t). In the linear case U(x) = x, Eq. (3) maps directly onto the Klein-Kramers FP ∂y 2∂y2 equation for the joint position-velocity PDF of a Brownian particle [1]. 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: Y˙(t)=F(y,t)+σ(y,t)ξ(t), (4) where ξ(t) is a Gaussian white noise with zero mean (cid:104)ξ(t)(cid:105) = 0 and covariance (cid:104)ξ(t)ξ(t(cid:48))(cid:105) = δ(t−t(cid:48)) and the Itˆo- convention is assumed for the multiplicative term σ(y,t)ξ(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) [1]. 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 ξ(t) with a L´evy noise, such that Y describes L´evy flight type dynamics, [21–25] 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 [26–32]. These fractional FK equations have been successfully used to model, e.g., the dynamics of migrating epithelial cells [33] and the advection of a fluid particle in turbulence [34]. However, the relation between such equations and the underlying stochastic dynamics is often not clear [35, 36]. Thus, more systematic approaches have been adopted, which explicitly assume the process Y to represent a continuous time randomwalk(CTRW)withjumplengthsandwaitingtimesdrawnfromindependentdistributions[37,38]. Specifically, in[39,40],startingfromarandomwalkdescriptionoftheCTRWinphasespace,aKlein-Kramersequationcontaining a fractional substantial derivative, which generalizes the ordinarymaterial derivative through the inclusionof explicit retardation effects, was derived. In [41] 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 [42, 43], as well as to inhomogeneous media in [44, 45]. 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 equa- tionsasintheconventionalpictureofEqs.(3,4). InsteadofusingaLangevin-typeequationtodescribethedynamics of the underlying CTRW, the time evolution of P(cid:98) is derived directly by means of a generalized master equation. Such full correspondence has been established only recently in our work [46] 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 allowsinparticulartocapturestraightforwardlydifferentwaitingtimedistributionsoftheCTRWbyamonotonically increasing L´evy process in an auxiliary time variable. The characteristic Laplace exponent Φ of the L´evy process is naturally related to the memory kernel appearing in the generalized FK equation. Our FK formula has been recently confirmed in [47], within a random walk approach, for the special case of tempered L´evy-stable distributed waiting times. As shown in [46], by employing a variable parametric form of Φ, one can fit the resulting anomalous process to mean-squaredisplacement(MSD)datadisplayinganonlinearcrossoverbetween,e.g.,subdiffusiveandnormaldiffusive 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 Φ, 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 [33, 48–51] to macromolecules and living organelles, e.g., mitochondria, inside the cytoplasm [52–63]. Thus, our framework can be applied to a large variety of different systems exhibiting such anomalous diffusive behavior. OurmainpurposeinthepresentmanuscriptistoextendthederivationofthegeneralizedFKequationforCTRWs withanarbitrarywaitingtimedistribution[46]tobothspace-andtime-dependentforcesandexplicittime-dependent functionals. These latter ones in particular, to our knowledge, have not yet been considered in previous works on 3 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 [43] within a random walk approach for power-law distributed waiting times, its extension to arbitrary distributions as expressed in the subordination approach of [46] 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´evy processes, semimartingales and their stochastic calculus, which is necessary to understandthemathematicalframeworkunderlyingthedescriptionofanomalousprocessesintermsofsubordination. 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 pedagogicalreasons, wefirsttreatthecaseofaspace-dependentforceandtime-independentfunctionalasin[46]. We thenextendthemethodtospace-andtime-dependentforcesandtime-dependentfunctionals. InSec.Vweapplyour 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 ThediscussionofrandomwalkswitharbitrarilydistributedwaitingtimesgoesbacktotheseminalworkbyMontroll and Weiss [37]. 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 [64–70]. A natural parametrisation of such a random walk is obtained in terms of the number n of jumps performed. If we call ξ the amplitude of the jump occurring at the jth step and by η the waiting time between the j j (j−1)th and jth jumps, the position Y and the elapsed time T are given by summing all such n RVs: n n (cid:88) (cid:88) Y =y + ξ , T = η , (5) 0 j j j=1 j=1 where y denotes the initial position. Rather than a parametrisation in terms of the discrete variable n, it is usually 0 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≥0:T(n)≤t}, and the position variable can be written as below: N(t) (cid:88) Y(t)=y + ξ . (6) 0 j j=1 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,t)=(cid:104)δ(y−Y(t))(cid:105) directly from Eq. (6). The master equation is then further approximated on a Y diffusive time and spatial scale leading to Fokker-Planck equations with fractional time derivatives, which describe the time evolution of f (y,t) [38, 71–76]. (ii) The diffusive limit can already be considered on the level of Eq. (5), Y thus leading to a coupled set of Langevin equations describing the stochastic process Y(t) [77–79]. The resulting Fokker-Planck equation for f (y,t) is equivalent to that obtained by using approach (i) [78, 80–82]. Y 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)→S(t) [83]. In such a continuum limit Eqs. (5) become: (cid:90) s (cid:90) s X(s)=y + ξ(s)ds, T(s)= η(s)ds, (7) 0 0 0 4 where s is now interpreted as an auxiliary or operational time variable. The position coordinate Eq. (6) becomes: (cid:90) S(t) Y(t)= ξ(τ)dτ =X(S(t)). (8) 0 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 S(t)= inf {s:T(s)>t}, (9) s>0 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. IIIC. 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: (cid:28)(cid:90) ∞ (cid:29) f (y,t)=(cid:104)δ(y−Y(t))(cid:105)= dsδ(s−S(t))δ(y−X(s)) Y 0 (cid:90) ∞ = dsh(s,t)f (y,s), (10) X 0 where f is the PDF associated with the process X(s) and h(s,t) that of S(t). X Regarding the waiting time process, a widely studied case is that of T being a one-sided L´evy-stable process of order 0 < α ≤ 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: (cid:68) (cid:69) e−λT(s) =e−sλα. (11) Thus, the time-change S in Eq. (9) is an inverse L´evy-stable subordinator with a PDF defined in Laplace space as (cid:101)h(s,λ)=(cid:82)∞dte−λth(s,t)=λα−1e−sλα [78]. When X describes a pure diffusion process with noise strength σ, the 0 MSD of Y exhibits a power-law scaling of the same exponent α, i.e., MSD(t)=(cid:10)[Y(t)−y ]2(cid:11)=σtα/Γ(1+α). This 0 particular scaling regime is characteristic of a pure subdiffusive system [38]. 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 [33, 48–57, 59–63]. 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´evy 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 [84, 85]: (cid:68) (cid:69) e−λT(s) =e−sΦ(λ), (12) with the Laplace exponent Φ characterizing the jump structure of the waiting times. Different functional forms of Φ correspond to different distribution laws of the waiting times and of the renewal process T. By choosing Φ suitably, severaldifferentwaitingtimestatisticscanbecaptured,i.e.,theanomalousprocessY(t)canbemodeledaccordingto the observed experimental dynamics. If we choose a power law Φ(λ)=λα, we recover Eq. (11), i.e., the CTRW case. If instead Φ(λ) = λ, T is a deterministic drift, T = s, and Y(t) reduces to a normal diffusion (Brownian limit) with exponentially distributed waiting times [38]. Details on the mathematical properties of Φ are discussed in Sec. IIIB. 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. (10). Indeed, the inverse of the process T has the PDF in Laplace space: (cid:101)h(s,λ)=[Φ(λ)/λ]e−sΦ(λ) [86]. AsthefirstandsecondmomentofX aregivenby(cid:104)X(s)(cid:105)=y0 and(cid:104)[X(s)]2(cid:105)= y02+σs respectively, those of the time-changed process Y read as (cid:104)Y(cid:101)(λ)(cid:105)=(cid:82)0∞(cid:101)h(s,λ)(cid:104)X(s)(cid:105)=y0/λ and (cid:104)Y(cid:101)2(λ)(cid:105)= (cid:82)∞(cid:101)h(s,λ)(cid:104)X2(s)(cid:105)=[y2+σ/Φ(λ)]/λ. Putting these results together, we obtain the MSD in Laplace space: 0 0 (cid:93) σ MSD(λ)= . (13) λΦ(λ) 5 We note that for Φ(λ) = λα we recover the single power-law scaling previously discussed. For different choices of Φ Eq. (13) is able to capture many different scaling behaviors of MSD(t) [46]. For instance, in the case of T given as a tempered L´evy-stable process, i.e., Φ(λ) = (µ+λ)α −µα with µ ∈ R+, the MSD displays a crossover between a power-law (of exponent α) and a normal linear scaling. In the case of a sum of two independent stable distributions with exponents α1 < α2, i.e., Φ(λ) = C1λα1 +C2λα2 with C1,C2 ∈ R+, the crossover is between two subdiffusive regimes with power-law scaling of exponents α /α respectively for long/small times [87]. 1 2 Equations (7) allow us to express X and T in terms of Langevin equations, as first formulated by Fogedby [77]: X˙(s)=ξ(s), T˙(s)=η(s). (14) withinitialconditionsX(0)=Y(0)=y andT(0)=0. Eventhoughthesestepsseemsomewhatsuperfluous,thereare 0 variousadvantagesbyexpressingbothX andT inthisway. Inparticular, itallowstoeasilyincorporateforcesacting on the random walker during the instantaneous jumps. The case where they affect the dynamics of the walker also duringthewaitingtimesisdiscussedin[88,89]. ThefactorizationoftheexpectedvaluesleadingtoEq.(10)stillholds inthepresenceofforcesthatdependonlyonthepositionorontheauxiliarytimevariable, i.e., whenEq.(14)(left)is substituted by X˙(s)=F(X(s),s)+ξ(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 [90]. Instead, the correct Langevin equation substituting Eq. (14)(left) is given as X˙(s)=F(X(s),T(s))+ξ(s) [79, 81, 82, 86, 90– 93]. Now the factorization in Eq. (10) 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. a) b) η ∆s 2 X(s) ∆x ∆x ∆x X(s) x∆ x∆ x ∆s x ∆ ∆ ∆s ∆s ∆s η3∆s x∆ ∆s t1 t2 t3 t4 T(s) η1∆s t1 t2 t3 t4 T(s) Figure1. IllustrationofstochasticprocessesdefinedbyasubordinationasinEqs.(8,9). Weconsider(i)adiscretisationofthe physical time t (black solid lines) of step length ∆t and (ii) a discretisation of finite step length ∆s (cid:28) ∆t of s (not shown). We denote with ∆x, ∆T the increments of the processes X and T corresponding to an increment ∆s (red dots), which are given by Eqs. (14). The resulting process Y(t) = X(S(t)), with S defined by Eq. (9), is plotted in black solid lines. The process X(s) is plotted in dashed red lines and considered as an ordinary random walk. ξ(s) is thus a white Gaussian noise. (a) Normal Diffusion. In this case T is a deterministic drift and ∆T = ∆s. Therefore, s coincides with the physical time t and its discretisation is a thinner time partition. (b) Anomalous Diffusion. In this case, ∆T = η ∆s, with η being a RV. j j Consequently, s no longer coincides with the physical time, but it provides a parametrisation of the elapsed physical time T viaEq.(9). Wenotetheoccurrenceoftrappingevents(Y(t )=Y(t )inpanelb),whichareduetothevariablelengthof∆T. 2 3 By specifying the properties of X and T, we can capture the fluctuation properties of a large variety of physi- cal systems. In the following, we generally assume ξ(s) to represent a white Gaussian noise with (cid:104)ξ(s)(cid:105) = 0 and (cid:104)ξ(s)ξ(s(cid:48))(cid:105) = δ(s−s(cid:48)), such that X is a normal diffusive process in the auxiliary time s. For generality, we also consider a multiplicative noise strength σ(x), which can capture, e.g., the effects of geometrical confinement (see [94] and references therein). Thus, we consider the following set of coupled Langevin equations: X˙(s)=F(X(s),T(s))+σ(X(s))ξ(s), T˙(s)=η(s), (15) wherethefunctionsF(x,t)andσ(x)satisfystandardconditions[95],weadopttheItˆoprescriptionforthemultiplica- tive 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). WenotethatinsteadofdefiningT directlybyEq.(12), itisconvenienttospecifytheunderlyingnoiseprocessη(s) in Eq. (15)(right) by means of its characteristic functional [84, 96]: G[k(s)]=(cid:68)e−(cid:82)0∞k(s)η(s)ds(cid:69)=e−(cid:82)0∞Φ(k(s))ds, (16) 6 which specifies the statistics of the whole noise trajectory. Recalling Eq. (7), we recover the characteristic function (cid:10)e−λT(s)(cid:11)fromEq.(16)bysettingk(s(cid:48))=λΘ(s−s(cid:48)). Thisresult,togetherwithEq.(15)(right),elucidatestherenewal (cid:82)∆s nature of the process T. Such a process is indeed expressed as a sum over waiting time increments ∆t= η(τ)dτ 0 over a small time step ∆s of characteristic function (cid:10)e−λ∆t(cid:11) = e−∆sΦ(λ), which can be used to simulate the process Y(t) within a suitable discretisation scheme [97]. Specifying the noise η 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´evy processes, semimartingales and time-changed processes. In particular, we characterize the function Φ, which specifies the waiting time statistics of the underlying random walk. III. MATHEMATICAL BASICS: LE´VY PROCESSES, SUBORDINATORS AND TIME-CHANGED PROCESSES In this section, we provide a comprehensive and accessible review of the theory of L´evy processes, subordinators andtime-changedprocesses,thatarenecessarytoformulateandworkwithLangevinequationsofanomalousdiffusive processes. For more details beyond our present discussion we refer to the monographs [84, 85]. A. L´evy processes A stochastic process Y(t) for t≥0 and initial condition Y(0)=y is a L´evy process if the conditions below hold: 0 1. Y(0)=y =0 almost surely (a.s.), i.e., for each of its different realizations. 0 2. Y(t) has independent increments, i.e., ∀n ≥ 2 and for each partition 0 ≤ t < t < ... < t ≤ t the RVs 0 1 n {Y(t )−Y(t )} are independent. j j−1 j=1,...,n 3. Y(t) has stationary increments, meaning that for all 0 ≤ t < t ≤ t the RV Y(t )−Y(t ) has the same 1 2 2 1 distribution as Y(t −t ). Note that, if 1 is not satisfied, it would instead depend on Y(t −t )−y . 2 1 2 1 0 4. The trajectories of Y(t) are c`adl`ag, i.e., right continuous with left limits. If one restricts the conditions 2, 4 by assuming Gaussian distributed increments and continuous trajectories respec- tively, one recovers ordinary Brownian motion. Moreover, as a consequence of (ii), Y is infinitely divisible ∀t≥0. The notion of infinitely divisibility characterizes a RV Y that can be expressed as a sum of different i.i.d. RVs {X }. Specifically, we assume Y to have PDF f and corresponding characteristic function: φ (k) = (cid:10)eikY(cid:11) = j Y Y (cid:82)+∞eikxf (x)dx. If ∀n ∈ N there exist i.i.d RVs X(n),...,X(n) with PDF f and characteristic function φ −∞ Y 1 n X X (uniquely defined), such that the following relation holds (in distribution): n Y =(cid:88)X(n), (17) j j=1 then Y is said to be infinitely divisible. Thus, their characteristic functions are related by the following equation: n φY(k)=(cid:68)eik(cid:80)nj=1Xj(n)(cid:69)= (cid:89)(cid:68)eikXj(n)(cid:69)=(cid:104)(cid:68)eikX1(n)(cid:69)(cid:105)n =[φX(k)]n. (18) j=1 Note that the factorization of the average is due to the independence of the RVs X(n) and that Eq. (18) represents j a necessary and sufficient condition for Y to be infinitely divisible [85], i.e., it can be used as a criterion to asses the infinitely divisibility of a given RV. Considering now the L´evy process Y in continuous time, we find that ∀n∈N and ∀t≥0 it can be rewritten as: n (cid:18) (cid:19) (cid:18) (cid:19) Y(t)=(cid:88)X(n), X(n) =Y jt −Y (j−1)t , (19) j j n n j=1 where {X(n)} are i.i.d. RVs because of 2-3. Therefore, such property uniquely relates the characteristic function of Y j to that of a general infinitely divisible RV. Indeed, let us introduce the auxiliary function: Ψ(k,t)=ln(cid:104)exp[ikY(t)](cid:105). (20) 7 By suitably adapting Eq. (19), we can write for arbitrary m,n∈N: m (cid:88) Y(m)=Y(1)+[Y(2)−Y(1)]+...+[Y(m)−Y(m−1)]= [Y(j)−Y(j−1)], (21a) j=1 n (cid:16)m(cid:17) (cid:104) (cid:16) m(cid:17) (cid:16)m(cid:17)(cid:105) (cid:104) (cid:16) m(cid:17)(cid:105) (cid:88)(cid:104) (cid:16) m(cid:17) (cid:16) m(cid:17)(cid:105) Y(m)=Y + Y 2 −Y +...+ Y(m)−Y (n−1) = Y j −Y (j−1) . (21b) n n n n n n j=1 Further exploiting the property 2, Eqs. (21a, 21b) imply that Y(m)=mY(1) and Y(m)=nY(cid:0)m(cid:1) (in distribution), n such that Ψ(k,m) can be computed exactly. Specifically, we obtain the two equivalent equations: Ψ(k,m)=ln(cid:104)exp[ikmY(1)](cid:105)=mln(cid:104)exp[ikY(1)](cid:105), (22a) (cid:68) (cid:104) (cid:16)m(cid:17)(cid:105)(cid:69) (cid:68) (cid:104) (cid:16)m(cid:17)(cid:105)(cid:69) Ψ(k,m)=ln exp iknY =n ln exp ikY . (22b) n n Here, the factorization of the ensemble average is allowed by the independence of the increments. Consequently, the equations Ψ(k,m)=mΨ(k,1) and Ψ(k,m)=nΨ(k,m/n) hold simultaneously. If we combine them, we obtain: (cid:16) m(cid:17) m Ψ k, = Ψ(k,1). (23) n n This relation expresses the characteristic function of Y at the finite time t=m/n in terms of its value at time t=1. Asitissatisfiedforeveryintegerm,n,italsoholdsforanyrealpositivet,i.e.,Ψ(k,t)=tΨ(k,1). Moreover,according to Eq. (19), Y(1) is an infinitely divisible RV. Therefore, a L´evy process Y(t) has the characteristic function: φ (k,t)=(cid:104)exp[ikY(t)](cid:105)=exp[G(k)t], (24) Y where we set G(k) = Ψ(k,1), with Ψ(k,1) being fully specified as the logarithm of the characteristic function of a generalinfinitelydivisibleRV.Suchquantity,alongwiththecharacteristicfunctionofageneralL´evyprocessaccording to Eq. (24), is uniquely characterized by means of the L´evy-Khintchine representation. Introducing parameters b∈R and σ ≥0, the L´evy-Khintchine representation states that the characteristic function of any infinitely divisible RV is of the form φ(k)=exp[G(k)] with the function G defined below [85]: 1 (cid:90) G(k)=ibk− σk2+ [eiky−1−iky1 (y)]Π(dy), (25) 2 |y|<1 R/{0} with 1 (y) = 1 for y ∈ A or 1 (y) = 0 otherwise. In Eq. (25) Π is a so called L´evy measure, i.e., a probability A A measure satisfying the following condition: (cid:90) Max(|y|2,1)Π(dy)<∞. (26) R/{0} Thus,anyL´evyprocessY isuniquelycharacterizedbythetriplet(b,σ,Π),whichdeterminesitscharacteristicfunction trough Eqs. (24, 25). Two examples of L´evy processes are of fundamental importance: • Brownian motion with drift. In this case, the characteristic function is given by Eq. (24) with the specification: σ2 G(k)=ibk− k2. (27) 2 As we can write G as below: (cid:20) (cid:18) b σ2 (cid:19)(cid:21)n G(k)= exp ik − k2 =[G (k)]n, (28) n 2n X where we define G (k) = exp[ikb/n−k2σ2/(2n)], we deduce that G(k) is the characteristic function of an X infinitely divisible RV. The RVs Y(n) in Eq. (17) are Gaussian distributed with mean b/n and variance σ2/n. j • Compound Poisson process. A Compound Poisson process Y(t) on the interval [0,t] is defined as N(t) (cid:88) Y(t)= ξ , (29) j j=1 8 where the ξ are i.i.d. RVs with law f and N(t) is a Poisson process characterized by an intensity λ, i.e., j ξ p(N(t) = n) = exp(−λt)[(λt)n/n!]. This represents a pure jump process where jumps of length ξ , drawn as j i.i.d. RVs, occur at time points that are spaced by an exponentially distributed waiting time. Within this picture, λ 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: (cid:42)  N(t) (cid:43) (cid:42)(cid:42)  n (cid:12)(cid:12) (cid:43)(cid:43) (cid:88) (cid:88) (cid:12) φ(k,t)= expik ξj = expik ξj(cid:12)N(t)=n (cid:12) j=1 j=1 (cid:12) (cid:42)(cid:20)(cid:90) +∞ (cid:21)n(cid:43) = dyeikyf (y) ξ −∞ (cid:88)∞ (cid:20)(cid:90) +∞ (cid:21)n (λt)n =exp(−λt) dyeikyf (y) ξ n! n=0 −∞ (cid:20) (cid:18)(cid:90) +∞ (cid:19)(cid:21) =exp λt dyeikyf (y)−1 . (30) ξ −∞ Therefore, Y(t) is a L´evy process with (cid:90) +∞ G(k)=λ (cid:0)eiky−1(cid:1) f (y)dy. (31) ξ −∞ Thus, G(k) is again the characteristic function of an infinitely divisible RV, as Eq. (17) is satisfied by choosing the RVs Y(n) to have the characteristic function Eq. (31) with λ→λ/n. j Comparing these two examples with the L´evy-Khintchine Eq. (25) we note that L´evy processes can be interpreted intuitively as consisting of three different contributions: A deterministic drift, a continuous normal diffusion and a discontinuousjumpprocess. TheCompoundPoissonprocessrepresentsthesimplestexampleofthejumpcontribution, where the L´evy measure of the jump amplitudes is just a normalizable PDF: Π(dy) = λf (y)dy. The specifications ξ giveninEqs.(25,26)extendthiscasetoawiderclassofjumpprocesses,whichhavepossiblynon-normalizablelength distribution or may possess an infinite intensity of small jumps [98]. For physical applications, important examples of such more exotic processes are L´evy-stable and tempered L´evy-stable processes [99]. AL´evy-stableprocessisaL´evyprocesswithstabledistributedincrements,i.e.,G(k)inEq.(24)isthecharacteristic 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 {Yj}j=1,...,n. If real-valued sequences of parameters {cn}n∈N and {dn}n∈N exist, such that the following relation holds in distribution: n (cid:88) Y =c Y +d , (32) j n n j=1 then Y is called a stable RV. If d =0, then Y is strictly stable. From this definition, it is straightforward to see that n (i) Y is infinitely divisible [simply set X(n) =(Y −d /n)/c in Eq. (17)] and that (ii) the existence of Y represents j j n n a generalization of the central limit theorem. Indeed, Eq. (32) equivalently states that the sequences of partial sums √ {Sn}n∈N with Sn =(Y1+...+Yn−dn)/cn converge in distribution to Y. With the choice cn =σ n and dn =nm, this is the ordinary central limit theorem and Y is Gaussian distributed with mean m and variance σ2. For different choicesofc andd , weobtaininsteadageneralizedcentrallimittheorem[100]. However, theonlypossiblechoiceto n n satisfy Eq. (32) is given by c =σn1/α, with 0<α≤2, also called index of stability of the stable distribution [101]. n As stable distributions are infinitely divisible, their characteristic function is completely determined by Eq. (25). In particular,wehavetwopossiblecharacteristics: (i)(b,σ2,0)forα=2,implyingthatY isGaussian(meanb,variance σ2) and (ii) (b,0,Π) for α(cid:54)=2 with Π specified by the following formula (for c ,c ≥0 and c +c >0): 1 2 1 2 (cid:26)c y−1−α dy y ∈[0,∞) Π(dy)= 1 (33) c |y|−1−α dy y ∈(−∞,0) 2 9 By suitably changing coordinates in Eq. (25) [102], we obtain the following characterization of G(k): 1 G(k)=iµk− σ2k2, α=2 (34a) 2 (cid:104) (cid:16)πα(cid:17)(cid:105) G(k)=iµk−σα|k|α 1−iβ sign(k)tan , α(cid:54)=1,2 (34b) 2 (cid:20) (cid:21) 2 G(k)=iµk−σ|k| 1+iβ sign(k)log(|k|) , α=1 (34c) π for µ∈R, σ ≥0 and −1≤β ≤1. If X is a symmetric stable RV (β =0), then the function G reduces to: G(k)=−ρα|k|α, 0<α≤2 (35) √ with ρ=σ for 0<α<2 and ρ=σ/ 2 for α=2. Of particular importance for us are one-sided monotonically increasing L´evy processes, which can be used to implement a random time change. These processes are called subordinators. B. Subordinators We define a subordinator a one-dimensional L´evy process that is a.s. non-decreasing. Thus, if T(t) for t ≥ 0 is a subordinator, the following properties hold a.s.: (i) T(t) ≥ 0, ∀t ≥ 0 and (ii) T(t ) ≤ T(t ), ∀t ≤ t . According to 1 2 1 2 the discussion in Sec. IIIA, its characteristic function is determined by Eqs. (24, 25) for a subclass of characteristic triplets (b,σ,ν) that we need to determine. We first note that, if X(t) is a Brownian motion of variance σ2, we have: p(X(t) ≥ 0) = 1/2 = p(X(t) ≤ 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´evy symbol, i.e., σ = 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≥0 and Π(−∞,0)=0. Taking these requirements into account, the characteristic function of a general subordinator T is given as [103]: (cid:68) (cid:69) e−λT(t) =e−tΦ(λ), (36) (cid:90) +∞ Φ(λ)=bλ+ (1−e−λy)Π(dy), (37) 0 (cid:82)+∞ where one needs to further assume that Max(y,1)Π(dy) < ∞. Φ is the Laplace exponent of the subordinator. 0 AssuggestedinSec.II,ΦdeterminesthecharacteristicfunctionalEq.(16)ofthenoiseη appearinginEq.(15)(right). Weremarkthatonlytwoparametersdefineitsform, i.e., thecharacteristicsofT aredeterminedbytheduplet(b,Π). Using Eq. (24) and Jensen’s inequality, one can show that Φ(λ) must be a continuous, non negative, non decreasing and concave function. We also remark that Φ(0) = 0. In general, one can prove that Φ is a Bernstein function [104, 105]. Specific examples of subordinators are reviewed in the following: • L´evy stable subordinator. A subordinator T is L´evy stable if it has characteristic duplet (0,Π) with α Π(dy)= y−1−αdy. (38) Γ(1−α) If we substitute it inside Eq. (37), we obtain the following Laplace exponent: α (cid:90) ∞ λ (cid:90) ∞ Φ(λ)= (1−e−λy)y−1−αdy = e−λyy−αdy =λα. (39) Γ(1−α) Γ(1−α) 0 0 • Tempered L´evy stable subordinator. A subordinator T is tempered L´evy stable if it has characteristic duplet (0,Π) with the L´evy measure [84]: α Π(dy)= e−µyy−1−αdy µ>0. (40) Γ(1−α) 10 If we substitute it inside Eq. (37), we obtain the following Laplace exponent: α (cid:90) ∞ Φ(λ)= (1−e−λy)e−µyy−1−αdy Γ(1−α) 0 α (cid:20) (cid:90) ∞ (cid:90) ∞ (cid:21) = − (1−e−µy)y−αdy+ (1−e−(λ+µ)y)y−αdy =(λ+µ)α−µα, (41) Γ(1−α) 0 0 wherewesolvedtheintegralsbyemployingEq.(39). TheresultingprocessT interpolatesbetweenaL´evystable (of order parameter α) 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 λ→0(∞) limit of its Laplace transform Eq. (36), orequivalentlyofthefunctionΦderivedinEq.(41). Specifically,forλ→∞wefind: Φ(λ)∼λα,whichrecovers thecaseofaL´evystablesubordinator[seeEq.(39)]. Forλ→0instead,weobtain: Φ(λ)∼αµα−1λ,suchthat we can approximate the characteristic function of T as (cid:104)e−λT(t)(cid:105) ∼ 1−αµα−1tλ ∼ (1+αµα−1tλ)−1, which is the Laplace transform of an exponential distribution. C. Semimartingales and the stochastic calculus of time-changed processes LetusconsideraL´evyprocessX(t)andasubordinatorT(t). Thankstoitsmonotonicity,T canbeemployeddirectly asarandomparametrisationoftimedefininganewtime-changedprocessY troughtherelation: Y(t)=X(T(t)). Such processiseasilyshowntostillbeaL´evyprocess[85]. However,asdiscussedinSec.II,thisisnotthesituationarising forCTRWs,wheretheirrepresentationbycoupledLangevinequationsinthediffusivelimitinvolvestheinverseofthe processT, i.e., theprocessS definedinEq.(9). ThecrucialpointisthatS isgenerallynot a L´evy process. Rather, it ispartofamoregeneralclassofprocessescalledsemimartingales, whichalsocontainL´evyprocessesasaspecialcase (see below). An important theorem tracing back to the work of Jacod [106] states that semimartingales (and thus L´evyprocesses)subordinatedbyproperlydefinedtime-changes, e.g., theprocessS, areagainsemimartingales. 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) [95, 107], 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 . The process M is a martingale if the following relation on its conditional average holds [98]: s (cid:104)M(t)|M(s)=m (cid:105)=m . (42) s s It is instead called a sub-martingale if (cid:104)M(t)|M(s)=m (cid:105) ≥ m or a super-martingale if (cid:104)M(t)|M(s)=m (cid:105) ≤ m . s s s 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: Y(t)=M(t)+A(t) (43) where M(t) and A(t) are a martingale and a finite variation process with c`adla`g paths respectively. We recall that stochastic integration with respect to semimartingales is well defined [108]. For the sake of our discussion, we will only present their Itˆo formula. In the specific case of Y being a continuous semimartingale, this is given by: (cid:90) t 1(cid:90) t f(Y(t))−f(Y )= f(cid:48)(Y(τ))dY(τ)+ f(cid:48)(cid:48)(Y(τ))d[Y,Y] , (44) 0 2 τ 0 0 where [Y,Y] is the quadratic variation of Y (see Appendix C for a review). The extension of Eq. (44) to a M- t dimensional semimartingale Z reads as: (cid:88)M (cid:90) t 1 (cid:88)M (cid:90) t f(Z(t))−f(Z )= f(cid:48)(Z(τ))dZ(i)(τ)+ f(cid:48)(cid:48) (Z(τ))d[Z(i),Z(j)] , (45) 0 i 2 i,j τ i=1 0 i,j=1 0 where[Z(i),Z(j)] isthejointquadraticvariationofZ(i),Z(j),whichisdefinedanalogouslytothequadraticvariation t by substituting the squared increment in Eq. (C1) with the product of the increments of the two processes. We note that both [Y,Y] and [Z(i),Z(j)] are continuous increasing processes. The joint one also has finite variation paths t t [108]. For further properties of semimartingales and their theory of stochastic integration we refer to [108]. 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