Universal spectral features of different classes of random-diffusivity processes

Diffusive processes came to the attention of the broader scientific community with the experiments on 'active molecules' by Brown, who reported the jittery motion of granules of '1/4000th to 1/5000th of an inch in length' contained in pollen grains as well as control experiments on powdered inorganic rocks [1]. In the mid-19th century physician-physiologist Fick published his studies on salt fluxes between reservoirs of different concentrations connected by tubes [2]. To quantify the observed dynamics Fick introduced the diffusion equation ('Fick's second law') for the spatio-temporal concentration profile. A major breakthrough was the theoretical description of 'Brownian motion' and the diffusion equation in terms of probabilistic arguments by Einstein [3], Smoluchowski [4], and Sutherland [5]. Concurrently Pearson introduced the notion of the 'random walk' [6], and Langevin proposed the intuitive picture of the random force and the stochastic Langevin equation [7].

More recently, major advances in experimental techniques such as superresolution microscopy continue to provide unprecedented insight into the motion of submicron and even fluorescently tagged molecular tracers in complex environments such as living biological cells [8–11]. Concurrently, simulations are becoming ever more powerful and reveal the molecular dynamics in systems such as lipid membranes [12] or internal protein motion [13]. The data resulting from such complex systems unveil a number of new phenomena in the stochastic particle motion and thus call for new theoretical concepts [14–16] on top of already known approaches [17–19].

Among these new insights is that endogenous and introduced tracers in living biological cells perform anomalous diffusion of the form $\left\langle {\mathbf{r}}^{2}\left(t\right)\right\rangle \simeq {K}_{\alpha }{t}^{\alpha }$ in a wide range of systems [8,9,20]. For instance, subdiffusion with 0 < α < 1 was measured for messenger RNA probes in bacteria cells [21,22], for DNA loci and telomeres in bacteria and eukaryotic cells [22–24], for granules in yeast and human cells [25,26], as well as for the stochastic motion of biological membrane constituents [27,28]. In these cases the slower than Brownian, passive tracer motion is effected by the highly crowded nature of the environment, as can be studied in in vitro systems [29,30]. In fact, even small green fluorescent proteins of some 2 nm in size were shown to subdiffuse [31]. Conversely, superdiffusion with 1 < α < 2 in biological cells is caused by active motion of molecular motors due to consumption of biochemical energy units. Examples include the motor motion itself [32,33], the transport of introduced plastic beads in fibroblast cells [34], RNA cargo in neuron cells [35], and of granules in amoeba [36].

However, even when the mean squared displacement seemingly suggests Brownian motion based on the observation that α = 1, remarkable effects have been reported recently. Thus, the motion of micron-sized tracer beads moving along nanotubes as well as in entangled polymer networks was shown to be 'Fickian' (α = 1) yet the measured displacement distribution exhibited significant deviations from the expected Gaussian law: namely, an exponential distribution of the form P(r, t) ∝ exp(−|r|/λ(t)) with λ(t) ∝ t^1/2 was observed [37,38]. Similar 'Fickian yet non-Gaussian' diffusion was found for the tracer dynamics in hard sphere colloidal suspensions [39], for the stochastic motion of nanoparticles in nanopost arrays [40], of colloidal nanoparticles adsorbed at fluid interfaces [41–43] and moving along membranes and inside colloidal suspension [44], and for the motion of nematodes [45]. Even more complicated non-Gaussian distributions of displacements were recently observed in Dictyostelium discoideum cells [46,47] and protein-crowded lipid bilayer membranes [48]. While in some experiments the non-Gaussian shape of P(r, t) is observed over the entire experimental window, others report clear crossover behaviours from a non-Gaussian shape at shorter time scales to an effective Gaussian behaviour at longer time scales, for instance, see [37,38].

A non-Gaussian probability density along with the scaling exponent α = 1 of the mean squared displacement can be achieved in the superstatistical approach, in which it is assumed that individual Gaussian densities are averaged over a distribution of diffusivities [49–53]. A microscopic realisation of such a behaviour was proposed for a model of diffusion during a polymerisation process [54]. However, in superstatistics (and in the related process called generalised grey Brownian motion [55–57]) the distribution is a constant of the motion and thus no crossover behaviour as mentioned above can be described. In order to include such a non-Gaussian to Gaussian crossover models were introduced in which the diffusion coefficient is considered as a stochastic process itself. In this diffusing-diffusivity picture, originally proposed by Chubynsky and Slater [58], the stochastic dynamics of the diffusivity is characterised by a well-defined correlation time above which the diffusivity becomes equilibrated. Concurrently to this equilibration the ensuing form of P(r, t) becomes effectively Gaussian. Random-diffusivity models have since then been developed and analysed further, and their application is mainly the diffusive dynamics in heterogeneous systems [59–70].

In fact, stochastic models based on random diffusivities are ubiquitous in financial mathematics for the modelling of stock price dynamics. They are commonly known as stochastic volatility models and many different examples have been analysed in order to identify a proper description for the volatility [71]. Among them, one can find diffusion-based models, where the volatility is described with continuous sample paths, as well as more complicated dynamics where, for instance, jumps are also allowed or where the volatility is defined as a function of separate stochastic processes [72]. Financial mathematics hosts a rich variaty of random-diffusivity models, motivated by various aspects of the observed financial market data. Here we present a range of additional, new random-diffusivity models in the context of time series analysis, extending the range of available models beyond the diffusing-diffusivity model developed for Fickian yet non-Gaussian diffusion processes. These may in turn be useful for financial mathematics. As both fields are quickly expanding and new facets are being continuously unveiled, we are confident that the different models introduced here and their detailed features offer the necessary flexibility to account for the new observations.

The central purpose of our study here is twofold. First we analyse several new classes of random-diffusivity models, divided into two groups, jump models and functionals of Brownian motion. For both groups we consider several concrete examples and derive analytic solutions for the probability density function (PDF) Π(x, t). The PDF turns out to delicately depend on the precise formulation of the model: the central part may be Gaussian or non-Gaussian, and the tails may be of Gaussian, exponential, log-normal, or even power-law shape. The second goal we pursue here are the spectral properties of random-diffusivity processes. Namely, while earlier studies of the random-diffusivity dynamics were mainly concerned with the PDF and the mean squared displacement encoded in the process we assume a different stance and derive the spectral properties of single particle trajectories with finite observation time, geared for the description of contemporary single particle tracking experiments. Such an analysis was worked out in detail for specific systems of normal and anomalous diffusion [73–78], and we here study the commonalities and differences emerging for random-diffusivity scenarios.

Traditionally, power spectral analyses are based on the textbook definition of the spectral density

$$\begin{equation}\mu \left(f\right)=\underset{T\to \infty }{\mathrm{lim}}\frac{1}{T}\left\langle {\left\vert {\int }_{0}^{T}{\mathrm{e}}^{\mathrm{i}ft}x\left(t\right)\mathrm{d}t\right\vert }^{2}\right\rangle .\end{equation} \tag{ 1 }$$

This definition involves taking the limit of infinite (practically, very long) measurement times as well as averaging over an ensemble (practically, a large number of) of particles, here and in the following denoted by angular brackets, $\left\langle \cdot \right\rangle$. Typical single particle tracking experiments, however, are limited in the measurement time, for instance, due to the lifetime of the employed fluorescent tags or the time a particle stays in the microscope focus. At the same time, such experiments are often limited to a relatively small number of individual trajectories. To cater for this common type of experimental situations we avoid taking the long time and ensemble limits by considering the single-trajectory power spectral density (PSD)

$$\begin{equation}{S}_{T}\left(f\right)=\frac{1}{T}{\left\vert {\int }_{0}^{T}{\mathrm{e}}^{\mathrm{i}ft}x\left(t\right)\mathrm{d}t\right\vert }^{2}\end{equation} \tag{ 2 }$$

as functions of frequency f and measurement time T. We previously analysed the behaviour of S_T(f) for different diffusion scenarios [76–78] and demonstrated that it is practically useful in the analysis of experimental data [76,77]. In what follows we derive the moment-generating function (MGF) of the PSD (2) for different classes of random-diffusivity processes, including several cases not yet studied in literature. In particular, we obtain the probability density P(A) of the single PSD amplitude, an intrinsically random quantity for a finite-time measurement of a stochastic motion that was demonstrated to be a very useful quantity for the analysis of measured particle trajectories. In addition to analytical derivations we present detailed numerical analyses. This study provides a quite general approach to obtain the PDF for any diffusing-diffusivity model, providing new insights on this class of processes.

This work is structured as follows. We start from section 2 with a description of the model and in section 3 we report general results on the spectral properties of this class of processes. Specific examples of diffusing-diffusivity models are described in sections 4–6. The first example is the well known case in which the diffusivity is modelled as the squared Ornstein–Uhlenbeck process. In the second group of examples we analyse two cases in which the diffusivity is defined as a jump process. The third and last group shows three examples in which the diffusivity is described as a functional of Brownian motion. Finally, in section 7 we draw our conclusions. In the appendix we report details on the explicit derivations of our results.

We consider a class of one-dimensional stochastic processes x_t that obey the Langevin equation in the Itô convention,

$$\begin{equation}{\dot {x}}_{t}=\sqrt{2{D}_{0}{{\Psi}}_{t}}{\xi }_{t}.\end{equation} \tag{ 3 }$$

Here D₀ is a constant, dimensional coefficient in units length²/time, and in our analysis we will assume the initial condition x₀ = 0. In equation (3) ξ_t denotes a standard Gaussian white noise with zero mean and covariance function $\overline{{\xi }_{t}{\xi }_{{t}^{\prime }}}=\delta \left(t-{t}^{\prime }\right)$. The bar here and henceforth denotes averaging with respect to the noise ξ_t. Lastly Ψ_t is a positive-definite random function, which multiplies D₀ and thus introduces a time-dependent randomness into the effective noise amplitude. In the following we stipulate that Ψ_t is Riemann-integrable on a finite interval (0, T) such that ${\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}$ exists with probability 1. Note that the case Ψ_t ≡ 1 corresponds to standard Brownian motion, while a deterministic choice of the form Ψ_t = t^α−1 produces so-called scaled Brownian motion [79,80]. We will here discuss several particular choices for the random function Ψ_t. In addition to the previously made choice of a squared Ornstein–Uhlenbeck process we will consider the case when Ψ_t is a jump-process, that attains independent, identically distributed random values. We also present several examples when Ψ_t is subordinated to standard unbiased Brownian motion B_t: namely, ${{\Psi}}_{t}={B}_{t}^{2}/{a}^{2}$, where a is a model parameter, Ψ_t = Θ(B_t), where Θ(x) is the Heaviside theta function, and geometric Brownian motion Ψ_t = exp(−B_t/a).

Regardless of the choice of the random function Ψ_t, we can solve the Langevin equation (3) for the trajectory x_t for a fixed realisation of the noise and a given realisation of Ψ_t, to obtain

$$\begin{equation}{x}_{t}={\left(2{D}_{0}\right)}^{1/2}{\int }_{0}^{t}\mathrm{d}\tau {{\Psi}}_{\tau }^{1/2}{\xi }_{\tau }.\end{equation} \tag{ 4 }$$

The characteristic function of x_t can be written down in the form

$$\begin{equation}{{\Phi}}_{w}={\left\langle \overline{\mathrm{exp}\left(\mathrm{i}w{\left(2{D}_{0}\right)}^{1/2}{\int }_{0}^{t}\mathrm{d}\tau {{\Psi}}_{\tau }^{1/2}{\xi }_{\tau }\right)}\right\rangle }_{{\Psi}},\end{equation} \tag{ 5 }$$

where the bar stands for averaging over thermal histories, while the angular brackets denote averaging over the realisations of the random function Ψ_t. The thermal average can be performed straightforwardly to give

$$\begin{equation}{{\Phi}}_{w}={\left\langle \mathrm{exp}\left(-{D}_{0}{w}^{2}{\int }_{0}^{t}\mathrm{d}\tau {{\Psi}}_{\tau }\right)\right\rangle }_{{\Psi}}.\end{equation} \tag{ 6 }$$

The desired PDF Π(x, t) can then be written as

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{d}w{\mathrm{e}}^{-\mathrm{i}wx}{{\Phi}}_{w}.\end{equation} \tag{ 7 }$$

In the following section 4 we provide several examples with explicit expressions for the probability density, and we will see how different choices of Ψ_t may lead to PDFs of considerably different shapes.

We first obtain exact expressions for the PSD (2) and then study the limiting behaviour for high frequencies.

3.1. Exact expressions for arbitrary frequency and observation time

We investigate the PSD of an individual trajectory x_t encoded in the stochastic dynamics (3) with t ∈ (0, T),

$$\begin{equation}{S}_{T}\left(f\right)=\frac{1}{T}{\int }_{0}^{T}\mathrm{d}{t}_{1}{\int }_{0}^{T}\mathrm{d}{t}_{2}\;\mathrm{cos}\left(f\left({t}_{1}-{t}_{2}\right)\right){x}_{{t}_{1}}{x}_{{t}_{2}},\end{equation} \tag{ 8 }$$

as function of the frequency f and the observation time T. We determine the MGF and the PDF of the random variable S_T(f).

The MGF of the single-trajectory PSD in (8) is defined as

$$\begin{equation}{\phi }_{\lambda }={\left\langle \overline{\mathrm{exp}\left(-\frac{\lambda }{T}{\int }_{0}^{T}\mathrm{d}t{\int }_{0}^{T}\mathrm{d}{t}^{\prime }\;\mathrm{cos}\left(f\left(t-{t}^{\prime }\right)\right){x}_{t}{x}_{{t}^{\prime }}\right)}\right\rangle }_{{\Psi}}\end{equation} \tag{ 9 }$$

with λ ⩾ 0. Relegating some intermediate calculations to appendix A we find the following expression for ϕ_λ in (9) averaged over thermal noises,

$$\begin{equation}{\phi }_{\lambda }=\frac{1}{4\pi \lambda }{\int }_{-\infty }^{\infty }\mathrm{d}{z}_{1}{\int }_{-\infty }^{\infty }\mathrm{d}{z}_{2}\;\mathrm{exp}\left(-\frac{{z}_{1}^{2}+{z}_{2}^{2}}{4\lambda }\right){\left\langle \mathrm{exp}\left(-{D}_{0}{\int }_{0}^{T}\;\mathrm{d}t{{\Psi}}_{t}{\left({\int }_{t}^{T}\mathrm{d}\tau {Q}_{\tau }\right)}^{2}\right)\right\rangle }_{{\Psi}},\end{equation} \tag{ 10 }$$

where

$$\begin{equation}{Q}_{t}={z}_{1}\frac{\mathrm{cos}\left(ft\right)}{\sqrt{T}}+{z}_{2}\frac{\mathrm{sin}\left(ft\right)}{\sqrt{T}}.\end{equation} \tag{ 11 }$$

Performing the inverse Laplace transform of expression (10) we find the general result for the PDF

$$\begin{equation}p\left({S}_{T}\left(f\right)=S\right)=\frac{1}{4\pi }{\int }_{-\infty }^{\infty }\mathrm{d}{z}_{1}{\int }_{-\infty }^{\infty }\mathrm{d}{z}_{2}{J}_{0}\left(\sqrt{\left({z}_{1}^{2}+{z}_{2}^{2}\right)S}\right){\left\langle \mathrm{exp}\left(-{D}_{0}{\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}{\left({\int }_{t}^{T}\mathrm{d}\tau {Q}_{\tau }\right)}^{2}\right)\right\rangle }_{{\Psi}},\end{equation} \tag{ 12 }$$

where J₀(z) denotes the Bessel function of the first kind. A more explicit dependence on the frequency f can be obtained in the form (see appendix A for more details)

$$\begin{align}\hfill {\phi }_{\lambda }& =\left\langle \left[1+\frac{8\lambda {D}_{0}}{{f}^{2}T}{\int }_{0}^{T}\;\mathrm{d}t{{\Psi}}_{t}\left(1-\mathrm{cos}\left(f\left(T-t\right)\right)\right)\right.\right.\hfill \\ \hfill & {\left.{\left.\quad +\;\frac{16{\lambda }^{2}{D}_{0}^{2}}{{f}^{4}{T}^{2}}{\int }_{0}^{T}\;\mathrm{d}{t}_{1}{{\Psi}}_{{t}_{1}}{\int }_{0}^{T}\;\mathrm{d}{t}_{2}{{\Psi}}_{{t}_{2}}\left(\frac{3}{4}+{L}_{f}\left({t}_{1},{t}_{2}\right)\right)\right]}^{-1/2}\right\rangle }_{{\Psi}},\hfill \end{align} \tag{ 13 }$$

where L_f(t₁, t₂) is defined by the somewhat lengthy expression (A.6).

The expression within the angular brackets in relation (13) is the exact MGF of the PSD of the process x_t in (3) for any fixed realisation of Ψ_t and holds for arbitrary T and arbitrary f. It also represents the exact form of the MGF in the case when Ψ_t is non-fluctuating: as mentioned, in particular, for Ψ_t = 1 it describes the MGF in case of standard Brownian motion [76], while the choice Ψ_t = t^α−1 corresponds to the case of scaled Brownian motion recently studied in [78].

3.2. Exact high frequency limiting behaviour

As already remarked we here concentrate on random processes Ψ_t which, for any finite T, are Riemann-integrable with probability 1, which implies that in the limit f → ∞ certain integrals vanish, as shown in appendix B. As a consequence, expression (13) attains the following exact analytic high-frequency form

$$\begin{equation}{\phi }_{\lambda }\sim {\left\langle {\left[1+\frac{8\lambda {D}_{0}}{{f}^{2}T}{\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}+\frac{12{\lambda }^{2}{D}_{0}^{2}}{{f}^{4}{T}^{2}}{\left({\int }_{0}^{T}\;\mathrm{d}t\;{{\Psi}}_{t}\right)}^{2}\right]}^{-1/2}\right\rangle }_{{\Psi}},\end{equation} \tag{ 14 }$$

in which we dropped the vanishing terms and kept only the leading terms in 1/f.

We note that the Laplace parameter λ appears in the combination D₀λ/f² so that the high-f spectrum of a single-trajectory PSD has the universal form

$$\begin{equation}{S}_{T}\left(f\right)\sim \frac{4{D}_{0}A}{{f}^{2}},\end{equation} \tag{ 15 }$$

regardless of the specific choice of Ψ_t. Here A is a dimensionless, random amplitude, which differs from realisation to realisation. This means that the characteristic high-frequency dependence of the PSD can be learned, in principle, from just a single trajectory, in agreement with the conclusions in [76–78].

The MGF Φ_λ of the random amplitude A follows from (14) and can be written as

$$\begin{align}\hfill {{\Phi}}_{\lambda }& ={\int }_{0}^{\infty }\mathrm{d}A{\mathrm{e}}^{-\lambda A}P\left(A\right)\hfill \\ \hfill & =\frac{2}{\sqrt{3}}{\int }_{0}^{\infty }\mathrm{d}p\;\mathrm{exp}\left(-\frac{4p}{3}\right){I}_{0}\left(\frac{2p}{3}\right){\Upsilon}\left(T;\lambda p/T\right),\hfill \end{align} \tag{ 16 }$$

where I₀(z) is the modified Bessel function of the first kind, and

$$\begin{equation}{\Upsilon}\left(T;\lambda \right)={\left\langle \mathrm{exp}\left(-\lambda {\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}\right)\right\rangle }_{{\Psi}}\end{equation} \tag{ 17 }$$

is the MGF of the integrated diffusivity (see [64,66])

$$\begin{equation}{\tau }_{T}={\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}.\end{equation} \tag{ 18 }$$

Relation (16) links the MGFs of A and τ_T. Moreover, note that the characteristic function of the diffusing-diffusivity model in (6) is tightly related to the MGF of τ_T in (17), specifically

$$\begin{equation}{\Upsilon}\left(T;{D}_{0}{w}^{2}\right)={{\Phi}}_{w}.\end{equation} \tag{ 19 }$$

As shown in [64] the function ϒ(T; λ) determines the first-passage time properties of the stochastic process x_t. Here we show how this function controls the high-frequency behaviour of the PSD.

Taking the inverse Laplace transform with respect to the parameter λ we evaluate the PDF of A,

$$\begin{equation}P\left(A\right)=\frac{2}{\sqrt{3}}{\int }_{0}^{\infty }\mathrm{d}z{J}_{0}\left(\left(1+\frac{1}{\sqrt{3}}\right)\sqrt{2zA}\right){J}_{0}\left(\left(1-\frac{1}{\sqrt{3}}\right)\sqrt{2zA}\right){\Upsilon}\left(T;z/T\right).\end{equation} \tag{ 20 }$$

This exact expression determines the high-f behaviour of the PDF p(S_T(f) = S),

$$\begin{equation}p\left({S}_{T}\left(f\right)=S\right)\sim \frac{{f}^{2}}{4{D}_{0}}P\left(A=\frac{S{f}^{2}}{4{D}_{0}}\right)\end{equation} \tag{ 21 }$$

as f → ∞. The exact high-f forms in (16) and (20) will serve as the basis of our analysis for several particular choices of the process Ψ_t in section 4.

Before proceeding, we stop to make several general statements.

• (a)  
  Expanding the exponential function on the right-hand side of (16) into the Taylor series in powers of λ we obtain straightforwardly the relation between the moments of A and the moments of the integrated diffusivity τ_T, which is valid for any n,

  $$\begin{equation}\mathbb{E}\left\{{A}^{n}\right\}={\left(\frac{3}{4}\right)}^{n+1/2}n!{}_{2}F_{1}\left(\frac{n+1}{2},\frac{n+2}{2};1;\frac{1}{4}\right){\left\langle \frac{{\tau }_{T}^{n}}{{T}^{n}}\right\rangle }_{{\Psi}},\end{equation} \tag{ 22 }$$

  where ${}_{2}F_{1}\left(a,b;c;z\right)$ is the Gauss hypergeometric function. Since the moments of τ_T are related to the moments of the process x_T [66] we also find

  $$\begin{equation}\mathbb{E}\left\{{A}^{n}\right\}={\left(\frac{3}{4}\right)}^{n+1/2}\frac{{\left(n!\right)}^{2}}{\left(2n\right)!}\;{}_{2}F_{1}\left(\frac{n+1}{2},\frac{n+2}{2};1;\frac{1}{4}\right){\left\langle \frac{{x}_{T}^{2n}}{{\left({D}_{0}T\right)}^{n}}\right\rangle }_{{\Psi}}.\end{equation} \tag{ 23 }$$
• (b)  
  Starting from the results in (22) and (23) we can readily obtain the moments of S_T(f) as well. In particular, if we focus on its average value, we have

  $$\begin{equation}\left\langle {S}_{T}\left(f\right)\right\rangle =\frac{4{D}_{0}{C}_{1}}{{f}^{2}T}{\left\langle {\tau }_{T}\right\rangle }_{{\Psi}}=\frac{2{C}_{1}}{{f}^{2}T}{\left\langle \overline{{x}_{T}^{2}}\right\rangle }_{{\Psi}},\end{equation} \tag{ 24 }$$

  where ${C}_{1}={\left(3/4\right)}^{3/2}{}_{2}F_{1}\left(1,3/2;1;1/4\right)$. This suggests that those random-diffusivity models that display anomalous scaling of the MSD, i.e., ${\left\langle \overline{{x}_{T}^{2}}\right\rangle }_{{\Psi}}\not\simeq T$, exhibit ageing behaviour, namely, a dependence of the PSD properties on the trajectory length T.
• (c)  
  Equations (15) and (22) permit us to directly access the coefficient of variation γ of the PDF p(S_T(f) = S) in the high-f limit. We get straightforwardly

  $$\begin{align}\hfill \gamma & ={\left(\frac{{\left\langle \overline{{S}_{T}^{2}\left(f\right)}\right\rangle }_{{\Psi}}-{\left\langle \overline{{S}_{T}\left(f\right)}\right\rangle }_{{\Psi}}^{2}}{{\left\langle \overline{{S}_{T}\left(f\right)}\right\rangle }_{{\Psi}}^{2}}\right)}^{1/2}\approx {\left(\frac{{\left\langle {A}^{2}\right\rangle }_{{\Psi}}-{\left\langle A\right\rangle }_{{\Psi}}^{2}}{{\left\langle A\right\rangle }_{{\Psi}}^{2}}\right)}^{1/2}\hfill \\ \hfill & ={\left(\frac{3}{4}\frac{{\left\langle {x}_{T}^{4}\right\rangle }_{{\Psi}}}{{\left\langle {x}_{T}^{2}\right\rangle }_{{\Psi}}^{2}}-1\right)}^{1/2}={\left(\frac{9}{4}\frac{{\left\langle {\tau }_{T}^{2}\right\rangle }_{{\Psi}}}{{\left\langle {\tau }_{T}\right\rangle }_{{\Psi}}^{2}}-1\right)}^{1/2}\hfill \\ \hfill & ={\left(\frac{9}{4}\frac{{\partial }_{\lambda }^{2}{\Upsilon}\left(T;\lambda \right){\vert }_{\lambda =0}}{{\left({\partial }_{\lambda }{\Upsilon}\left(T;\lambda \right){\vert }_{\lambda =0}\right)}^{2}}-1\right)}^{1/2},\hfill \end{align} \tag{ 25 }$$

  which implies that the effective broadness of p(S_T(f) = S) is entirely defined by the first two moments of the random variable τ_T in (18). Specifically, it is independent of D₀ and f when f is only large enough.
• (d)  
  The behaviour of the left tail of p(S_T(f) = S) can be assessed in the following way. Note that the product of the two Bessel functions in (20) can be represented as a power series with an infinite radius of convergence (see (C.1) in appendix C). Inserting the expansion in (C.1) in (20) and integrating over z we find

  $$\begin{equation}P\left(A\right)=\frac{2}{\sqrt{3}}\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^{n}}{n!}{\left(\frac{\sqrt{3}+1}{\sqrt{6}}\right)}^{2n}{}_{2}F_{1}\left(-n,-n;1;\frac{1-\sqrt{3}/2}{1+\sqrt{3}/2}\right){\left\langle \frac{{T}^{n+1}}{{\tau }_{T}^{n+1}}\right\rangle }_{{\Psi}}{A}^{n},\end{equation} \tag{ 26 }$$

  if the inverse moments of the variable τ_T exist (and do not grow too fast with n). Therefore, the PDF P(A) is an analytic function of A in the vicinity of A = 0, with

  $$\begin{equation}P\left(0\right)=\frac{2}{\sqrt{3}}{\left\langle \frac{T}{{\tau }_{T}}\right\rangle }_{{\Psi}}.\end{equation} \tag{ 27 }$$

  We note that below we will encounter both situations when P(A) is analytic and when it is not. In the latter situation we will show that p(S_T(f) = S) diverges as S → 0, which can be already inferred from (27).

In this and the following sections we apply the above general theory to several random-diffusivity models. According to our main results (16) and (20) one first needs to evaluate the MGF ϒ(T; λ) of the integrated diffusivity τ_T for a chosen diffusivity process Ψ_t. To illustrate the quality of the theoretical predictions in the high-frequency limit we also performed numerical simulations using a Python code. The Euler integration scheme is used to compute (3), where Ψ_t is obtained by a numerical integration of the proper stochastic equation for each case. The PSD is obtained by fast Fourier transform for each trajectory. Starting from the single-trajectory power spectra the random amplitude A is calculated according to (15).

Concretely when Ψ_t in the diffusing-diffusivity model is defined as a stochastic process satisfying some Langevin equation, the distribution of A is determined by (16) and (20) through the MGF ϒ(T; λ) of the integrated diffusivity τ_T that can be obtained by solving the associated backward Fokker–Planck equation (see [66] for details). Here we consider the common example of squared Ornstein–Uhlenbeck process and related models.

The Ornstein–Uhlenbeck process Y_t defined by the stochastic equation

$$\begin{equation}{\dot {Y}}_{t}=-{\tau }_{\star }^{-1}{Y}_{t}+{\sigma }_{\star }{\xi }_{t}^{\prime }\end{equation} \tag{ 28 }$$

is a stationary Gaussian process mean-reverting to zero at a time scale τ_⋆ and driven by standard Gaussian white noise ξ'_t with volatility σ_⋆. The process ${{\Psi}}_{t}={Y}_{t}^{2}$ is one of the most common models of diffusing-diffusivity, which satisfies, due to Itô's formula,

$$\begin{equation}{\dot {{\Psi}}}_{t}={\tau }^{-1}\left(\overline{{\Psi}}-{{\Psi}}_{t}\right)+\sigma \sqrt{2{{\Psi}}_{t}}{\xi }_{t}^{\prime },\end{equation} \tag{ 29 }$$

where τ = τ_⋆/2, $\sigma =\sqrt{2}{\sigma }_{\star }$, and $\overline{{\Psi}}={\sigma }_{\star }^{2}{\tau }_{\star }/2={\sigma }^{2}\tau /2$. This model was extended in [59–61] by considering Ψ_t as the sum of n independent squared Ornstein–Uhlenbeck processes, when (29) still holds with $\overline{{\Psi}}=n{\sigma }^{2}\tau /2$. More generally, setting $\overline{{\Psi}}$ to be any positive constant, the Langevin equation (29) defines the so-called Feller process [81], also known as square root process or the Cox–Ingersoll–Ross process [82], and related in the Heston model [83]. This process was used to model the diffusing-diffusivity in [63,64], see also the discussion in [61].

The MGF ϒ(T; λ) for the integrated squared Ornstein–Uhlenbeck process was first computed by Dankel [84] and employed in [59–61]. Its computation for the Feller process in (29) was presented in [63],

$$\begin{equation}{\Upsilon}\left(T;\lambda \right)={\left(\frac{4\omega {\mathrm{e}}^{-\left(\omega -1\right)T/\left(2\tau \right)}}{{\left(\omega +1\right)}^{2}-{\left(\omega -1\right)}^{2}{\mathrm{e}}^{-\omega T/\tau }}\right)}^{\nu },\end{equation} \tag{ 30 }$$

where $\omega =\sqrt{1+4{\sigma }^{2}{\tau }^{2}\lambda }$ and $\nu =\overline{{\Psi}}/\left(\tau {\sigma }^{2}\right)$. In particular, setting $\overline{{\Psi}}={\sigma }^{2}\tau /2$ (and thus ν = 1/2) one retrieves the MGF for the squared Ornstein–Uhlenbeck process. A detailed discussion on the PDF of this model is presented in [59–61,63,64]. Using the explicit formulas for ${\left\langle \overline{{x}_{T}^{2}}\right\rangle }_{{\Psi}}$ and ${\left\langle \overline{{x}_{T}^{4}}\right\rangle }_{{\Psi}}$ from [61,63] we get from (25) that

$$\begin{equation}\gamma ={\left[\frac{3}{4}\left(3+\frac{6\tau }{\nu T}\left(1-\frac{\tau }{T}\left(1-{\mathrm{e}}^{-T/\tau }\right)\right)\right)-1\right]}^{1/2}.\end{equation} \tag{ 31 }$$

Moreover, as the second moment $\left\langle \overline{{x}_{T}^{2}}\right\rangle$ shows a linear trend in time [61,63], no ageing of the PSD occurs, as suggested in (24).

The PDF of A is determined via (20). Since an explicit calculation of this integral is not straightforward we perform a numerical integration. The results are shown in figure 1, in which we observe excellent agreement between the simulations and the theoretical results. The 1/f² scaling is recovered for any value of τ_⋆. The coefficient of variation γ converges to different values when we change τ_⋆, according to (31). Note that this result reflects the different degrees of broadness of the PDF of the random amplitude A. In particular for τ_⋆ ≪ T we obtain a result that is very similar to the one of Brownian motion, while for increasing τ_⋆ the PDF of the random amplitude A becomes increasingly broader.

Zoom In Zoom Out Reset image size

We divide the interval (0, T) into N equal subintervals of duration δ = T/N and suppose that Ψ_t is a jump process on these intervals, of the form

$$\begin{equation}{{\Psi}}_{t}={\psi }_{k}\;\;\text{on}\;\;t\in \left(\left[k-1\right]\delta ,k\delta \right),\quad k=1,\dots ,\;N.\end{equation} \tag{ 32 }$$

Furthermore we stipulate that the ψ_k are independent, identically distributed, positive-definite random variables with PDF ρ(ψ). In other words, we take that Ψ_t at each discrete time instant (k − 1)δ attains a new random value, taken from the common distribution, and stays constant and equal to this value up to the next discrete instant kδ. For a given realisation of the process Ψ_t we thus have

$$\begin{equation}{\tau }_{T}={\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}=\delta \sum _{k=1}^{N}{\psi }_{k},\end{equation} \tag{ 33 }$$

and hence

$$\begin{equation}{\Upsilon}\left(T;\lambda \right)={\left({\int }_{0}^{\infty }\mathrm{d}\psi \rho \left(\psi \right){\mathrm{e}}^{-\lambda \delta \psi }\right)}^{T/\delta }.\end{equation} \tag{ 34 }$$

Evaluating explicitly the derivatives ∂_λϒ(T; λ) and ${\partial }_{\lambda }^{2}{\Upsilon}\left(T;\lambda \right)$ at λ = 0, we get

$$\begin{equation}\left\langle {\tau }_{T}\right\rangle =T\;\mathbb{E}\left\{{\psi }_{k}\right\},\end{equation} \tag{ 35 }$$

when the first moment $\mathbb{E}\left\{{\psi }_{k}\right\}$ exists. From this we infer ${\left\langle \overline{{x}_{T}^{2}}\right\rangle }_{{\Psi}}$ and thus the respective ageing behaviour. Moreover, the coefficient of variation becomes

$$\begin{equation}\gamma ={\left[\frac{9}{4}\left(1-\frac{\delta }{T}+\frac{\delta }{T}\frac{\mathbb{E}\left\{{\psi }_{k}^{2}\right\}}{\mathbb{E}{\left\{{\psi }_{k}\right\}}^{2}}\right)-1\right]}^{1/2},\end{equation} \tag{ 36 }$$

when the first two moments $\mathbb{E}\left\{{\psi }_{k}\right\}$ and $\mathbb{E}\left\{{\psi }_{k}^{2}\right\}$ exist.

Modelling the diffusivity as a jump process can be seen as a way to describe the model in section 4 through a different parametrisation. Indeed, we define a time scale, which is given by the duration δ of each step interval, and we then introduce a random variability of the diffusivity from one interval to the next. These diffusivity fluctuations are chosen according to the PDF ρ(ψ). Of course, the main difference comes from the fact that in this model we do not have any correlation between successive diffusivities. In what follows we analyse two examples in detail. In the first one we select a Gamma distribution for ρ(ψ), in analogy with the diffusing-diffusivity model in section 4, where the diffusion coefficient shows a Gamma distribution as well. In the second example we select a Lévy–Smirnov distribution for ρ(ψ). This allows us to model a system in which a high probability of having small values of the diffusivity is combined with the presence of few outliers, which can be related, for instance, to values of the diffusivity at boundaries of the system.

5.1. Example I: Gamma distribution

First, we consider the Gamma distribution,

$$\begin{equation}\rho \left(\psi \right)=\frac{{\psi }^{\nu -1}}{{\Gamma}\left(\nu \right){\psi }_{0}^{\nu }}\;\mathrm{exp}\left(-\psi /{\psi }_{0}\right)\end{equation} \tag{ 37 }$$

with the shape parameter ν > 0 and the scale parameter ψ₀ > 0. From (34), we deduce

$$\begin{equation}{\Upsilon}\left(T;\lambda \right)={\left(1+\lambda \delta {\psi }_{0}\right)}^{-\nu T/\delta }.\end{equation} \tag{ 38 }$$

Position-PDF Π(x, t)

A direct calculation of the PDF for this model can be performed. Starting from (7) and recalling that Φ_w = ϒ(T; D₀w²), we get

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)={\mathcal{N}}_{t}{\left(\frac{\vert x\vert }{2\sqrt{{D}_{0}\delta {\psi }_{0}}}\right)}^{\frac{\nu t}{\delta }-\frac{1}{2}}{K}_{\frac{1}{2}-\frac{\nu t}{\delta }}\left(\frac{\vert x\vert }{\sqrt{{D}_{0}\delta {\psi }_{0}}}\right),\end{equation} \tag{ 39 }$$

where the normalisation coefficient is

$$\begin{equation}{\mathcal{N}}_{t}=\sqrt{\frac{2}{\pi }}/\left({\left({D}_{0}\delta {\psi }_{0}\right)}^{3/4}{\Gamma}\left(\frac{\nu t}{\delta }\right)\right).\end{equation} \tag{ 40 }$$

With the properties

$$\begin{equation}{z}^{\nu }{K}_{-\nu }\left(z\right)\sim {2}^{\nu -1}{\Gamma}\left(\nu \right)-\frac{{2}^{\nu -3}{\Gamma}\left(\nu \right)}{\nu -1}{z}^{2}\end{equation} \tag{ 41 }$$

for |z| → 0 and ν > 1, as well as

$$\begin{equation}{K}_{-\nu }\sim \sqrt{\frac{\pi }{2z}}{\mathrm{e}}^{-z}\end{equation} \tag{ 42 }$$

for |z| → ∞, the asymptotic behaviours of the PDF are given by

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)\sim {\mathcal{N}}_{t}\;{2}^{\frac{\nu t}{\delta }-\frac{3}{2}}{\Gamma}\left(\frac{\nu t}{\delta }-\frac{1}{2}\right)\left[1-\frac{{x}^{2}}{4\left(\frac{\nu t}{\delta }-\frac{3}{2}\right){D}_{0}\delta {\psi }_{0}}\right]\end{equation} \tag{ 43 }$$

for |x| → 0 and νt > 3δ/2, as well as

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)\sim \frac{2}{\sqrt{{D}_{0}\delta {\psi }_{0}}{\Gamma}\left(\frac{\nu t}{\delta }\right)}{\left(\frac{\vert x\vert }{2\sqrt{{D}_{0}\delta {\psi }_{0}}}\right)}^{\frac{\nu t}{\delta }-1}\mathrm{exp}\left(-\frac{\vert x\vert }{\sqrt{{D}_{0}\delta {\psi }_{0}}}\right)\end{equation} \tag{ 44 }$$

for |x| → ∞.

The functional behaviour of the PDF Π(x, t) is shown in figure 2. We see that by changing δ we can observe different shapes of Π(x, t). When δ = 1 [panels (a) and (c)] the Gaussian approximation (51) already provides a good estimate of the PDF over a wide range. We start observing discrepancies only far out in the tails, for values which can hardly be reached with real data. When δ = 100 [panels (b) and (d)], in contrast, the exponential tails are distinct. The behaviours at small and large x are well described by the asymptotic expansions in (43) and (44). Note that the value of δ in here plays a role similar to the correlation time τ_⋆ in the diffusing-diffusivity model defined in section 4. The only difference is that by changing τ_⋆ in the model above we also change the average diffusivity while, in this case, changes in the value of δ do not affect the average diffusivity, which is fixed once we choose the jumps PDF in (37).

Zoom In Zoom Out Reset image size

Amplitude-PDF P(A)

The MGF of the amplitude A of the jump process-diffusivity model is given by

$$\begin{equation}{{\Phi}}_{\lambda }=\frac{2}{\sqrt{3}}{\int }_{0}^{\infty }\mathrm{d}p\frac{\mathrm{exp}\left(-4p/3\right){I}_{0}\left(2p/3\right)}{{\left(1+p{\psi }_{0}\lambda \delta /T\right)}^{\nu T/\delta }},\end{equation} \tag{ 45 }$$

so that

$$\begin{equation}P\left(A\right)=\frac{2}{\sqrt{3}}{\int }_{0}^{\infty }\mathrm{d}z\frac{{J}_{0}\left(\left(1+1/\sqrt{3}\right)\sqrt{2zA}\right)}{{\left(1+z{\psi }_{0}\delta /T\right)}^{\nu T/\delta }}{J}_{0}\left(\left(1-1/\sqrt{3}\right)\sqrt{2zA}\right).\end{equation} \tag{ 46 }$$

In particular one has $\mathbb{E}\left\{{\psi }_{k}\right\}=\nu {{\Psi}}_{0}$ and $\mathbb{E}\left\{{\psi }_{k}^{2}\right\}={{\Psi}}_{0}^{2}\nu \left(\nu +1\right)$, thus from (35) we readily obtain ${\left\langle \overline{{x}_{T}^{2}}\right\rangle }_{{\Psi}}\simeq T$, demonstrating that in this process no ageing behaviour is displayed. Moreover, from (36) we get

$$\begin{equation}\gamma ={\left[\frac{9}{4}\left(1+\frac{\delta }{\nu T}\right)-1\right]}^{1/2}.\end{equation} \tag{ 47 }$$

In the limit δ → 0 and N → ∞, with δN = T fixed, we have

$$\begin{equation}{\Upsilon}\left(T;\lambda \right)\sim \mathrm{exp}\left(-\nu {\psi }_{0}T\lambda \right).\end{equation} \tag{ 48 }$$

Hence,

$$\begin{align}\hfill {{\Phi}}_{\lambda }& =\frac{2}{\sqrt{3}}{\int }_{0}^{\infty }\mathrm{d}p\;\mathrm{exp}\left(-\left(\frac{4}{3}+\nu {\psi }_{0}\lambda \right)p\right){I}_{0}\left(\frac{2p}{3}\right)\hfill \\ \hfill & ={\left(1+2\nu {\psi }_{0}\lambda +\frac{4}{3}{\left(\nu {\psi }_{0}\lambda \right)}^{2}\right)}^{-1/2}\hfill \end{align} \tag{ 49 }$$

and

$$\begin{align}\hfill P\left(A\right)& =\frac{2}{\sqrt{3}}{\int }_{0}^{\infty }\mathrm{d}z{J}_{0}\left(\left(1+1/\sqrt{3}\right)\sqrt{2zA}\right){J}_{0}\left(\left(1-1/\sqrt{3}\right)\sqrt{2zA}\right){\mathrm{e}}^{-\nu {\psi }_{0}z}\hfill \\ \hfill & =\frac{2}{\sqrt{3}\nu {\psi }_{0}}\;\mathrm{exp}\left(-\frac{4A}{3\nu {\psi }_{0}}\right){I}_{0}\left(\frac{2A}{3\nu {\psi }_{0}}\right).\hfill \end{align} \tag{ 50 }$$

This means that we have essentially the same behaviour as for standard one-dimensional Brownian motion, however, with renormalised coefficients (compare with the result in [76]), in agreement also with what we obtained for the diffusing-diffusivity model in section 4, when selecting τ_⋆ ≪ T. Indeed, if we use (48) and recall that Φ_w = ϒ(T; D₀w²), we readily obtain

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)\sim \frac{1}{2\sqrt{\pi \nu {\psi }_{0}{D}_{0}t}}\;\mathrm{exp}\left(-\frac{{x}^{2}}{4\nu {\psi }_{0}{D}_{0}t}\right).\end{equation} \tag{ 51 }$$

In figure 3 we show a direct comparison between the numerical and theoretical results for the Gamma distribution with ψ₀ = 1 and ν = 0.5. We observe that the average value of the power spectrum is not affected by the value of δ. Nevertheless, when we plot some sample single-trajectory power spectra we notice a larger amplitude scatter for larger values of δ. This may be clearly seen in the distribution of the random variable A, which is broader for larger values of δ, and consequently in the different limiting values of the coefficient of variation. Thus, the fluctuations are sensitive to different parameters of the distribution (37), while the mean behaviour is not.

Zoom In Zoom Out Reset image size

5.2. Example II: Lévy-Smirnov distribution

In our second example we consider the Lévy–Smirnov distribution

$$\begin{equation}\rho \left(\psi \right)=\sqrt{\frac{{\psi }_{0}}{\pi }}\frac{\mathrm{exp}\left(-{\psi }_{0}/\psi \right)}{{\psi }^{3/2}},\end{equation} \tag{ 52 }$$

for which equation (34) yields

$$\begin{equation}{\Upsilon}\left(T;\lambda \right)=\mathrm{exp}\left(-2T\sqrt{{\psi }_{0}\lambda /\delta }\right).\end{equation} \tag{ 53 }$$

Note that in this case $\mathbb{E}\left\{{\psi }_{k}\right\}$ and $\mathbb{E}\left\{{\psi }_{k}^{2}\right\}$ are not defined, such that ${\left\langle \overline{{x}_{T}^{2}}\right\rangle }_{{\Psi}}$ does not exist either. This suggests that a clear ageing behaviour cannot be defined and that fluctuations are what dominates the system.

Position-PDF Π(x, t)

As a consequence, we obtain the following analytical expression for the PDF,

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)=\frac{2t\sqrt{{D}_{0}{\psi }_{0}/\delta }}{\pi }\frac{1}{4{t}^{2}{D}_{0}{\psi }_{0}/\delta +{x}^{2}},\end{equation} \tag{ 54 }$$

where we recognise the power-law behaviour, that is already built into relation (52). Note that expression (54) represents the Cauchy distribution, whose median grows with time t.

The PDF is shown for two different values of δ in figure 4. We observe that, differently from the case with the Gamma distribution above, we do not see significant changes in the shape of the distribution when varying δ. For both cases, δ = 1 and δ = 100, the power-law behaviour (54) is readily discernible.

Zoom In Zoom Out Reset image size

Amplitude-PDF P(A)

The MGF for the random amplitude A reads

$$\begin{equation}{{\Phi}}_{\lambda }=\frac{2}{\sqrt{3}}{\int }_{0}^{\infty }\mathrm{d}p\;\mathrm{exp}\left(-\frac{4}{3}p-2\sqrt{p{\psi }_{0}\lambda T/\delta }\right){I}_{0}\left(\frac{2p}{3}\right)\end{equation} \tag{ 55 }$$

and

$$\begin{align}\hfill P\left(A\right)& =\frac{2}{\sqrt{3}}{\int }_{0}^{\infty }\mathrm{d}z{J}_{0}\left(\left(1+1/\sqrt{3}\right)\sqrt{2zA}\right){J}_{0}\left(\left(1-1/\sqrt{3}\right)\sqrt{2zA}\right)\mathrm{exp}\left(-2\sqrt{z{\psi }_{0}T/\delta }\right)\hfill \\ \hfill & =\frac{\delta }{\sqrt{3}{\psi }_{0}T}\frac{1}{{\left(1+\xi \right)}^{3/2}}\;{}_{2}F_{1}\left(\frac{3}{4},\frac{5}{4};1;\frac{{\xi }^{2}}{4{\left(1+\xi \right)}^{2}}\right),\hfill \end{align} \tag{ 56 }$$

with ξ = (4Aδ)/(3ψ₀T). Note that in the limit A → ∞, the leading behaviour of P(A) follows

$$\begin{equation}P\left(A\right)\sim \frac{1}{{A}^{3/2}}.\end{equation} \tag{ 57 }$$

Thus, the PDF P(A) inherits the property of diverging moments from the parental Lévy–Smirnov distribution.

Figure 5 summarises the properties of the PSD for the jump process with Lévy–Smirnov distribution (ψ₀ = 0.5). We observe that, despite the fat-tailed PDF in (54) we still observe the universal 1/f² scaling of the PDF. Concurrently, the PDF of the random amplitude A features the power-law behaviour according to (56). Note that the non-existence of the moments of P(A) generates a pronounced scatter in the amplitude of the average power spectrum.

Zoom In Zoom Out Reset image size

We now focus on the case when Ψ_t is a genuine 'diffusing-diffusivity' in the sense that it is subordinated to Brownian motion B_t starting at the origin at t = 0, with zero mean and covariance function

$$\begin{equation}\left\langle {B}_{t}{B}_{{t}^{\prime }}\right\rangle =2{D}_{B}\mathrm{min}\left\{t,{t}^{\prime }\right\}.\end{equation} \tag{ 58 }$$

We here choose Ψ_t = V[B_t], where V is some prescribed, positive-definite function. Note that random variables of the form ${\int }_{0}^{T}\mathrm{d}tV\left[{B}_{t}\right]$ appear across many disciplines, including probability theory, statistical analysis, computer science, mathematical finance and physics. Starting from earlier works [87–92], much effort has been invested in the analysis of the PDF and the corresponding Laplace transforms of these processes. A large body of exact results has been obtained within the last seven decades (see, e.g., [93–98] and further references therein). In the following, we consider three particular examples of V[B_t], for which we can carry out exact calculations and obtain insightful results.

6.1. Example I: Ψ_t = Θ(B_t)

First, we choose the cut-off Brownian motion

$$\begin{equation}{{\Psi}}_{t}={\Theta}\left({B}_{t}\right),\end{equation} \tag{ 59 }$$

where Θ(x) is the Heaviside theta function. The process x_t exhibits standard diffusive motion, once B_t > 0, and pauses, remaining at the position it has reached when B_t goes to negative values. The random variable ${\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}$ defines the time spent by a Brownian trajectory, starting at the origin, on the positive real line within the time interval (0, T). The time intervals between any two 'diffusion tours', as well as their duration, are random variables with a broad distribution.

This example is of particular interest as it represents an alternative to other standard models describing waiting times and trapping events. One could think, for instance, of the comb model, where a particle, while performing standard Brownian motion along one direction, gets stuck for a random time in branches perpendicular to the direction of the diffusive motion [99,100].

The MSD of the process x_t, as one can straightforwardly check, is just

$$\begin{equation}{\left\langle \overline{{x}_{t}^{2}}\right\rangle }_{{\Psi}}={D}_{0}t,\end{equation} \tag{ 60 }$$

that is, a standard diffusion law in which the diffusion coefficient is reduced by the factor 1/2. This means that no ageing behaviour is observed. For higher order moments one expects, of course, significant departures from the standard diffusive behaviour.

Position-PDF Π(x, t)

The MGF of the random variable ${\tau }_{T}={\int }_{0}^{T}\mathrm{d}t{\Theta}\left({B}_{t}\right)$, which has a bounded support on (0, T), was first derived by Kac [90], and Erdös and Kac [91]. Rewriting their result in our notation we have

$$\begin{align}\hfill {\Upsilon}\left(T;\lambda \right)& ={\left\langle \mathrm{exp}\left(-\lambda {\int }_{0}^{T}\mathrm{d}t{\Theta}\left({B}_{t}\right)\right)\right\rangle }_{{\Psi}}\hfill \\ \hfill & ={\mathrm{e}}^{-\lambda T/2}{I}_{0}\left(\frac{\lambda T}{2}\right).\hfill \end{align} \tag{ 61 }$$

Note that the inverse Laplace transform of this expression produces the celebrated Lévy arcsine law [101]. With result (61), the desired PDF is given by

$$\begin{align}\hfill \mathbf{\Pi }\left(x,t\right)& =\frac{1}{\pi }{\int }_{0}^{\infty }\mathrm{d}w\;\mathrm{cos}\left(wx\right)\;\mathrm{exp}\left(-\frac{{D}_{0}t}{2}{w}^{2}\right){I}_{0}\left(\frac{{D}_{0}t}{2}{w}^{2}\right)\hfill \\ \hfill & =\frac{\mathrm{exp}\left(-{x}^{2}/\left(8{D}_{0}t\right)\right)}{2\sqrt{{\pi }^{3}{D}_{0}t}}{K}_{0}\left(\frac{{x}^{2}}{8{D}_{0}t}\right).\hfill \end{align} \tag{ 62 }$$

Recalling that K₀(z) ∼ −ln(z/2) − γ_EM for |z| → 0, where γ_EM is the Euler–Mascheroni constant, for small x we have

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)\sim \frac{-\mathrm{ln}\left(\frac{{x}^{2}}{8{D}_{0}t}\right)-{\gamma }_{\text{EM}}}{2\sqrt{{\pi }^{3}{D}_{0}t}},\quad \vert x\vert \to 0.\end{equation} \tag{ 63 }$$

For large x we use (42) and obtain the asymptotic behaviour for the PDF,

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)\sim \frac{1}{\pi \vert x\vert }\mathrm{exp}\left(-\frac{{x}^{2}}{4{D}_{0}t}\right),\quad \vert x\vert \to \infty .\end{equation} \tag{ 64 }$$

The PDF for this process is shown in figure 6. We see that the central part of the PDF is strongly non-Gaussian, while the tails are Gaussian, in agreement with the asymptotic behaviours (63) and (64).

Zoom In Zoom Out Reset image size

Amplitude-PDF P(A)

Inserting expression (61) into (16) and performing the integration over z, we arrive at the following, remarkably compact expression for the MGF,

$$\begin{equation}{{\Phi}}_{\lambda }=\frac{2\sqrt{2}}{\pi \sqrt{2+3\lambda }}\mathbf{K}\left(\frac{2\lambda }{2+3\lambda }\right)\end{equation} \tag{ 65 }$$

where K(x) is the complete elliptic integral of the first kind,

$$\begin{equation}\mathbf{K}\left(x\right)={\int }_{0}^{\pi /2}\frac{\mathrm{d}\phi }{\sqrt{1-x\;{\mathrm{sin}}^{2}\left(\phi \right)}}.\end{equation} \tag{ 66 }$$

Note that the high-f asymptotic form in (65) is independent of the observation time T.

To proceed we take advantage of the definition of the complete elliptic integral and perform the inverse Laplace transform of (65). After some formal manipulations this yields the following expression for the PDF,

$$\begin{equation}P\left(A\right)=2\sqrt{\frac{2}{3{\pi }^{3}A}}{\int }_{0}^{\pi /2}\frac{\mathrm{d}\phi }{\sqrt{1-\frac{2}{3}\;{\mathrm{sin}}^{2}\left(\phi \right)}}\;\mathrm{exp}\left(-\frac{2A}{3\left(1-\frac{2}{3}\;{\mathrm{sin}}^{2}\left(\phi \right)\right)}\right).\end{equation} \tag{ 67 }$$

Multiplying both sides of the latter equation by Aⁿ and integrating over A from 0 to ∞, we get the following simple expression for the moments of the random amplitude A of arbitrary order,

$$\begin{equation}\mathbb{E}\left\{{A}^{n}\right\}=\frac{{\Gamma}\left(n+\frac{1}{2}\right)}{\sqrt{\pi }}{\left(\frac{3}{2}\right)}^{n}\;{}_{2}F_{1}\left(-n,\frac{1}{2};1;\frac{2}{3}\right).\end{equation} \tag{ 68 }$$

Then, from (25), we readily get the coefficient of variation, $\gamma =\sqrt{19/8}$.

Note that the integrals ${\int }_{0}^{\infty }\mathrm{d}\lambda {\lambda }^{n}{\Upsilon}\left(T;\lambda \right)$ diverge for any n > 0, which means that τ_T does not have negative moments. One therefore expects that P(A) is a non-analytic, diverging function in the limit A → 0. The small-A asymptotic behaviour of P(A) can be deduced directly from (67). Expanding the exponential function in the integral into the Taylor series in powers of A and expressing the emerging generalised elliptic integrals via their representations in terms of the toroidal functions ${P}_{n-1/2}\left(\mathrm{cosh}\left(\eta \right)\right)$ (see (C.3) in appendix C), we get

$$\begin{equation}P\left(A\right)=\sqrt{\frac{2}{\sqrt{3}\pi A}}\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^{n}}{n!}{\left(\frac{2}{\sqrt{3}}\right)}^{n}{P}_{n-1/2}\left(\frac{2}{\sqrt{3}}\right){A}^{n}.\end{equation} \tag{ 69 }$$

For the opposite limit A → ∞, we conveniently rewrite equation (56) in the form

$$\begin{align}\hfill P\left(A\right)& =2\sqrt{\frac{2}{3{\pi }^{3}A}}\;\mathrm{exp}\left(-\frac{2A}{3}\right){\int }_{0}^{\pi /2}\frac{\mathrm{d}\phi }{\sqrt{\left(1-\frac{2}{3}\;{\mathrm{sin}}^{2}\left(\phi \right)\right)}}\;\mathrm{exp}\left(-\frac{2A}{3}\frac{\frac{2}{3}\;{\mathrm{sin}}^{2}\left(\phi \right)}{1-\frac{2}{3}\;{\mathrm{sin}}^{2}\left(\phi \right)}\right)\hfill \\ \hfill & =2\sqrt{\frac{2}{3{\pi }^{3}A}}\;\mathrm{exp}\left(-\frac{2A}{3}\right)\sum _{n=0}^{\infty }\left({\int }_{0}^{\pi /2}\;{\mathrm{sin}}^{2n}\left(\phi \right)\right){\left(\frac{2}{3}\right)}^{n}{\mathrm{L}}_{n}^{\left(-1/2\right)}\left(\frac{2A}{3}\right)\hfill \\ \hfill & =\frac{2}{3\pi }\sqrt{\frac{3}{2A}}\;\mathrm{exp}\left(-\frac{2A}{3}\right)\sum _{n=0}^{\infty }\frac{{\Gamma}\left(n+1/2\right)}{n!}{\left(\frac{2}{3}\right)}^{n}{\mathrm{L}}_{n}^{\left(-1/2\right)}\left(\frac{2A}{3}\right),\hfill \end{align} \tag{ 70 }$$

where ${\mathrm{L}}_{n}^{\left(-1/2\right)}\left(x\right)$ are associated Laguerre polynomials. We focus next on the asymptotic behaviour of the function

$$\begin{equation}g\left(u\right)={u}^{-1/2}\sum _{n=0}^{\infty }\frac{{\Gamma}\left(n+1/2\right)}{n!}{\left(\frac{2}{3}\right)}^{n}{\mathrm{L}}_{n}^{\left(-1/2\right)}\left(u\right)\end{equation} \tag{ 71 }$$

in the limit u → ∞. Performing a Laplace transform of g(u) we readily get

$$\begin{align}\hfill {\mathcal{L}}_{s}\left\{g\left(u\right)\right\}& ={\int }_{0}^{\infty }\mathrm{d}u\;\mathrm{exp}\left(-su\right)g\left(u\right)\hfill \\ \hfill & =\frac{1}{{s}^{1/2}}\sum _{n=0}^{\infty }\frac{{{\Gamma}}^{2}\left(n+1/2\right)}{{\left(n!\right)}^{2}}{\left(\frac{2\left(s-1\right)}{3s}\right)}^{n}\hfill \\ \hfill & =\frac{2}{{s}^{1/2}}\mathbf{K}\left(\frac{2\left(s-1\right)}{3s}\right),\hfill \end{align} \tag{ 72 }$$

where K is the complete elliptic integral defined in equation (66). In the limit s → 0 (corresponding to A → ∞),

$$\begin{equation}{\mathcal{L}}_{s}\left\{g\left(u\right)\right\}\sim \frac{2}{{s}^{1/2}}\mathbf{K}\left(-\frac{2}{3s}\right)=2{\int }_{0}^{\pi /2}\frac{\mathrm{d}\phi }{\sqrt{s+\frac{2}{3}\;{\mathrm{sin}}^{2}\left(\phi \right)}}.\end{equation} \tag{ 73 }$$

Inverting the Laplace transform and integrating over ϕ, we find

$$\begin{equation}g\left(u\right)\sim \frac{\pi }{\sqrt{u}}\;\mathrm{exp}\left(-\frac{u}{3}\right){I}_{0}\left(\frac{u}{3}\right)\to \sqrt{\frac{3}{2}}\frac{1}{u}.\end{equation} \tag{ 74 }$$

Thus, in the limit A → ∞, the leading behaviour of the PDF P(A) yields in the form

$$\begin{equation}P\left(A\right)\sim \frac{1}{\pi A}\sqrt{\frac{3}{2}}\;\mathrm{exp}\left(-\frac{2A}{3}\right).\end{equation} \tag{ 75 }$$

Figure 7 summarises the numerical results for this case. Again we observe excellent agreement with the theoretical results.

Zoom In Zoom Out Reset image size

6.2. Example II: Ψ_t = exp(−B_t/a)

As the second example, we link the process Ψ_t to so-called geometric Brownian motion. Random variables of this form have been widely studied in the mathematical finance literature (see, e.g., [93]). Within the latter domain, they emerge very naturally as representation of the solution of the celebrated Black–Scholes equation. Their time-averaged counterpart is related to the so-called asian options [102–104] (see also [105,106]) and also appears in different contexts in the analysis of transport phenomena in disordered media (see, e.g., [107–113]) as well as characterises some features of the melting transition of heteropolymers [114].

In our notation, we set

$$\begin{equation}{{\Psi}}_{t}=\mathrm{exp}\left(-\frac{{B}_{t}}{a}\right),\end{equation} \tag{ 76 }$$

where a is a parameter of unit length. In this case, x_t exhibits an anomalously strong superdiffusion such that

$$\begin{align}\hfill \left\langle \overline{{x}_{t}^{2}}\right\rangle & =2{D}_{0}{\int }_{0}^{t}\mathrm{d}\tau {\left\langle {{\Psi}}_{\tau }\right\rangle }_{{\Psi}}=2{D}_{0}{\int }_{0}^{t}\mathrm{d}\tau \;\mathrm{exp}\left({D}_{\mathrm{B}}\tau /{a}^{2}\right)\hfill \\ \hfill & =\frac{2{D}_{0}{a}^{2}}{{D}_{\mathrm{B}}}\left(\mathrm{exp}\left({D}_{\mathrm{B}}t/{a}^{2}\right)-1\right).\hfill \end{align} \tag{ 77 }$$

Note that when D_Bt/a² ≪ 1 we have $\left\langle \overline{{x}_{t}^{2}}\right\rangle \sim 2{D}_{0}t$. These results for the MSD demonstrate that for this model, in general, we observe ageing behaviour, though the latter may be hidden while analysing very short trajectories.

Position-PDF Π(x, t)

The Laplace transform of the time-averaged geometric Brownian motion in our notation reads [108–110]

$$\begin{align}\hfill {\Upsilon}\left(T;\lambda \right)& ={\left\langle \mathrm{exp}\left(-\lambda {\int }_{0}^{T}\;\mathrm{d}t\;\mathrm{exp}\left(-\frac{{B}_{t}}{a}\right)\right)\right\rangle }_{{\Psi}}\hfill \\ \hfill & =\frac{2a}{\sqrt{\pi {D}_{\mathrm{B}}T}}{\int }_{0}^{\infty }\mathrm{d}x\;\mathrm{exp}\left(-\frac{{a}^{2}{x}^{2}}{{D}_{\mathrm{B}}T}\right)\;\mathrm{cos}\left(2a\sqrt{\frac{\lambda }{{D}_{\mathrm{B}}}}\;\mathrm{sinh}\left(x\right)\right)\hfill \\ \hfill & =\frac{2}{\pi }{\int }_{0}^{\infty }\mathrm{d}x\;\mathrm{exp}\left(-\frac{{D}_{\mathrm{B}}T}{4{a}^{2}}{x}^{2}\right)\;\mathrm{cosh}\left(\frac{\pi x}{2}\right){K}_{ix}\left(2a\sqrt{\frac{\lambda }{{D}_{\mathrm{B}}}}\right),\hfill \end{align} \tag{ 78 }$$

where K_ix is the modified Bessel function of the second kind with purely imaginary index. As a consequence, the PDF is given by

$$\begin{align}\hfill \mathbf{\Pi }\left(x,t\right)& ={\left(\frac{2}{\pi }\right)}^{3/2}{\int }_{0}^{\infty }\mathrm{d}z\;\mathrm{exp}\left(-\frac{{D}_{B}t}{4{a}^{2}}{z}^{2}\right)\;\mathrm{cosh}\left(\frac{\pi z}{2}\right){\int }_{0}^{\infty }\mathrm{d}w\;\mathrm{cos}\left(wx\right){K}_{iz}\left(2a\vert w\vert \sqrt{\frac{{D}_{0}}{{D}_{\mathrm{B}}}}\right)\hfill \\ \hfill & =\frac{1}{2\sqrt{\pi {b}_{2}t\left({b}_{1}^{2}+{x}^{2}\right)}}\;\mathrm{exp}\left(-\frac{{\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}}^{2}\left(x/{b}_{1}\right)}{4{b}_{2}t}\right),\hfill \end{align} \tag{ 79 }$$

where

$$\begin{equation}{b}_{1}=2a\sqrt{\frac{{D}_{0}}{{D}_{\mathrm{B}}}},\quad {b}_{2}=\frac{{D}_{\mathrm{B}}}{4{a}^{2}}.\end{equation} \tag{ 80 }$$

Recalling that arcsinh(z) ∼ z for z → 0 and arcsinh(z) ∼ ln(2z) for z → ∞, we express the asymptotic behaviour of the PDF as

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)\sim \frac{1}{2\sqrt{\pi {b}_{2}t\left({b}_{1}^{2}+{x}^{2}\right)}}\;\mathrm{exp}\left(-\frac{{x}^{2}}{4{b}_{2}t{b}_{1}^{2}}\right)\end{equation} \tag{ 81 }$$

for |x| → 0 and

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)\sim \frac{1}{2\sqrt{\pi {b}_{2}t\left({b}_{1}^{2}+{x}^{2}\right)}}\;\mathrm{exp}\left(-\frac{{\mathrm{ln}}^{2}\left(2x/{b}_{1}\right)}{4{b}_{2}t}\right)\end{equation} \tag{ 82 }$$

for |x| → ∞.

The PDF is shown in figure 8. In this case, according to the asymptotic expansions (81) and (82) we observe that the central part of the PDF is approximately Gaussian, while the tails follow a log-normal shape.

Zoom In Zoom Out Reset image size

Amplitude-PDF P(A)

Evaluating explicitly ∂_λϒ(T; λ) and ${\partial }_{\lambda }^{2}{\Upsilon}\left(T;\lambda \right)$ at λ = 0, we get

$$\begin{equation}\gamma ={\left[\frac{3}{8}\left(3+{\mathrm{e}}^{{D}_{\mathrm{B}}T/{a}^{2}}\left(2+{\mathrm{e}}^{{D}_{\mathrm{B}}T/{a}^{2}}\right)\right)-1\right]}^{1/2}.\end{equation} \tag{ 83 }$$

Figure 9 summarises our numerical results, which show excellent agreement with the theoretical prediction.

Zoom In Zoom Out Reset image size

6.3. Example III: ${{\Psi}}_{t}={B}_{t}^{2}/{a}^{2}$

Finally, we consider squared Brownian motion

$$\begin{equation}{{\Psi}}_{t}=\frac{{B}_{t}^{2}}{{a}^{2}},\end{equation} \tag{ 84 }$$

where a is a parameter of unit length. Note that here the process x_t in (3) is superdffusive, such that its mean-squared displacement obeys

$$\begin{equation}{\left\langle \overline{{x}_{t}^{2}}\right\rangle }_{{\Psi}}=\frac{2{D}_{0}{D}_{\mathrm{B}}}{{a}^{2}}{t}^{2},\end{equation} \tag{ 85 }$$

and thus a pronounced ageing behaviour occurs.

This example, similarly to the diffusing-diffusivity model in section 4, defines the diffusivity as the squared of an auxiliary variable, though in this case the variable follows a Brownian motion instead of the OU process. This choice implies that there is no crossover time, in contrast to the standard diffusing-diffusivity model, and thus we obtain a model which is always non-stationary. In particular, we introduce a larger separation between small and large values of the diffusivity, which may be interpreted as non-linear effects of the heterogeneity. Note that, if we were to define a random duration δ of the intervals, this model could be linked to a correlated CTRW [85,86].

Position-PDF Π(x, t)

The Laplace transform of the PDF of integrated squared Brownian motion was first calculated in the classical paper by Cameron and Martin [88,89] (see also [90]). In our notation,

$$\begin{align}\hfill {\Upsilon}\left(T;\lambda \right)& ={\left\langle \mathrm{exp}\left(-\frac{\lambda }{{a}^{2}}{\int }_{0}^{T}\mathrm{d}t{B}_{t}^{2}\right)\right\rangle }_{{\Psi}}\hfill \\ \hfill & =\frac{1}{\sqrt{\mathrm{cosh}\left(\sqrt{4{D}_{\mathrm{B}}{T}^{2}\lambda /{a}^{2}}\right)}},\hfill \end{align} \tag{ 86 }$$

and the PDF Π(x, t) for this process is given by

$$\begin{align}\hfill \mathbf{\Pi }\left(x,t\right)& =\frac{1}{\pi }{\int }_{0}^{\infty }\frac{\mathrm{d}w\;\mathrm{cos}\left(wx\right)}{\sqrt{\mathrm{cosh}\left(w\sqrt{4{D}_{\mathrm{B}}{D}_{0}{t}^{2}/{a}^{2}}\right)}}\hfill \\ \hfill & =\frac{1}{\sqrt{2}ct}\;\text{sech}\left(\frac{\pi x}{ct}\right){P}_{\frac{ix}{ct}-\frac{1}{2}}\left(0\right),\hfill \end{align} \tag{ 87 }$$

where $c=2\sqrt{{D}_{0}{D}_{\mathrm{B}}}/a$ and P_ν(z) is the Legendre function of the first kind. The latter admits the representation

$$\begin{equation}{P}_{ix/ct-1/2}\left(0\right)=\sqrt{\pi }/\left[{\Gamma}\left(\frac{ix}{2ct}+\frac{3}{4}\right){\Gamma}\left(-\frac{ix}{2ct}+\frac{3}{4}\right)\right],\end{equation} \tag{ 88 }$$

such that

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)=\frac{\sqrt{\pi }}{ct\sqrt{2}}\frac{\text{sech}\left(\frac{\pi x}{ct}\right)}{{\Gamma}\left(\frac{ix}{2ct}+\frac{3}{4}\right){\Gamma}\left(-\frac{ix}{2ct}+\frac{3}{4}\right)}.\end{equation} \tag{ 89 }$$

In the limits,

$$\begin{equation}{\Gamma}\left(\frac{ix}{2ct}+\frac{3}{4}\right){\Gamma}\left(-\frac{ix}{2ct}+\frac{3}{4}\right)\sim {\Gamma}{\left(\frac{3}{4}\right)}^{2}\end{equation} \tag{ 90 }$$

for |x| → 0, and

$$\begin{equation}{\Gamma}\left(\frac{ix}{2ct}+\frac{3}{4}\right){\Gamma}\left(-\frac{ix}{2ct}+\frac{3}{4}\right)\sim \pi \sqrt{\frac{2\vert x\vert }{ct}}\;\mathrm{exp}\left(-\frac{\pi \vert x\vert }{2ct}\right)\end{equation} \tag{ 91 }$$

for |x| → ∞. As a consequence, the behaviour of the PDF for small x is approximately Gaussian, ≃exp(−const.x²), where const. is a constant that can be expressed via the polylogarithm function. Conversely,

$$\begin{equation}\mathbf{\Pi }\left(x,t\right)\sim \frac{1}{\pi ct\vert x\vert }\;\mathrm{exp}\left(-\frac{\pi \vert x\vert }{2ct}\right)\end{equation} \tag{ 92 }$$

for |x| → ∞. The shape of the PDF is shown in figure 10. We clearly observe that the central part is approximately Gaussian while the tails have an exponential trend, as expressed explicitly by the asymptote (92).

Zoom In Zoom Out Reset image size

Amplitude-PDF P(A)

Using (86) we then find that the MGF of the random amplitude A and the corresponding PDF are given by the integrals

$$\begin{equation}{{\Phi}}_{\lambda }\left(A\right)=\frac{2}{\sqrt{3}}{\int }_{0}^{\infty }\mathrm{d}p\frac{\mathrm{exp}\left(-4p/3\right){I}_{0}\left(2p/3\right)}{\sqrt{\mathrm{cosh}\left(\sqrt{4{D}_{\mathrm{B}}T\lambda p/{a}^{2}}\right)}}\end{equation} \tag{ 93 }$$

and

$$\begin{equation}P\left(A\right)=\frac{2}{\sqrt{3}}{\int }_{0}^{\infty }\mathrm{d}z\frac{{J}_{0}\left(\left(1+1/\sqrt{3}\right)\sqrt{2Az}\right)}{\sqrt{\mathrm{cosh}\left(\sqrt{4{D}_{\mathrm{B}}Tz/{a}^{2}}\right)}}{J}_{0}\left(\left(1-1/\sqrt{3}\right)\sqrt{2Az}\right).\end{equation} \tag{ 94 }$$

The series representation of P(A) in (94) can be found directly by taking advantage of expansion (C.2) in appendix C and our result (56), to yield

$$\begin{equation}P\left(A\right)=\frac{1}{2\sqrt{6\pi }}\frac{{a}^{2}}{{D}_{\mathrm{B}}T}\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^{n}{\Gamma}\left(n+1/2\right)}{{\left(n+1/4\right)}^{2}n!{\left(1+{\xi }_{n}\right)}^{3/2}}\;{}_{2}F_{1}\left(\frac{3}{4},\frac{5}{4};1;\frac{{\xi }_{n}^{2}}{4{\left(1+{\xi }_{n}\right)}^{2}}\right)\end{equation} \tag{ 95 }$$

with

$$\begin{equation}{\xi }_{n}=\frac{{a}^{2}A}{3{\left(n+1/4\right)}^{2}{D}_{\mathrm{B}}T}.\end{equation} \tag{ 96 }$$

Note that the integrals

$$\begin{equation}{\int }_{0}^{\infty }\mathrm{d}\lambda {\lambda }^{n}{\Upsilon}\left(T;\lambda \right)\end{equation} \tag{ 97 }$$

exist for any n > 0 and, hence, all negative moments of ${\int }_{0}^{T}\mathrm{d}t{B}_{t}^{2}$ exist, as well. As a consequence, P(A) is an analytic function of A. We immediately obtain the coefficient of variation from (25) as $\gamma =\sqrt{17}/2$. Figure 11 displays numerical results for which we observe excellent agreement with the theoretical expressions.

Zoom In Zoom Out Reset image size

Quite typically stochastic time series are evaluated in terms of the ensemble averaged MSD ${\left\langle \overline{{x}_{t}^{2}}\right\rangle }_{{\Psi}}$. It has the advantage that fluctuations are reduced due to the averaging over many individual trajectories. However, this is not always possible. Namely, for the by-now routine results from single particle tracking experiments typically rather few, finite time series are obtained. These are then evaluated in terms of the time-averaged MSD. While this quantity may also be averaged over the available individual trajectories, it is increasingly realised that the amplitude fluctuations between individual trajectories in fact harbours important quantitative information characteristic for a given stochastic process [14,15,115–117].

Similar to the consideration of time averaged MSDs for trajectories of finite measurement time T, we here analysed the single-trajectory PSD of stochastic trajectories characterised by random diffusivities. Following our previous work on standard Brownian motion [76], as well as fractional [77] and scaled Brownian motion [78], we here investigated the detailed behaviour of single-trajectory PSDs of a broad class of diffusing-diffusivity models. These have recently gained considerable attention as simple models for diffusion processes in heterogeneous media. We described a general procedure to obtain the position PDF Π(x, t) for all such models. The main ingredient in the calculation is the MGF of the integrated diffusivity, showing explicitly that different choices of the underlying diffusivity Ψ_t lead to distinctly different emerging behaviours, as summarised in table 1.

Table 1. Collection of the main results for the different random-diffusivity models with respective equation numbers. Next to the definitions of the models we refer to the MGF of the integrated diffusivity τ_T, the coefficient of variation γ, which can play the role of an indicator for each model, and the ageing behaviour.

Random-diffusivity model		MGF of τT	γ	Ageing
${{\Psi}}_{t}={Y}_{t}^{2}$, Yt = OU process		equation (30)	equation (31)	×
Ψt = ψk	ρ(ψ) gamma distr.	equation (38)	equation (47)	×
	ρ(ψ) Lévy–Smirnov distr.	equation (53)	not defined	not defined
Ψt = V[Bt]	V[Z] = θ(Z)	equation (61)	$\sqrt{19/8}$	×
	V[Z] = exp(−Z)	equation (78)	equation (83)	✓
	V[Z] = Z2	equation (86)	$\sqrt{17}$/2	✓

We started our discussion from the by-now well-established and widely studied model of 'diffusing-diffusivity', namely the case when Ψ_t is chosen as the squared of the Ornstein–Uhlenbeck process. We then discussed the second case, in which Ψ_t is defined as a jump process. The properties of this model depend strongly on the exact distribution chosen for the increment variables. We considered two examples, a Gamma distribution and a Lévy–Smirnov distribution. Finally, three cases in which the diffusivity is modelled as a functional of Brownian motion were discussed.

The main result of this work is that, regardless of the different properties of all these diffusing-diffusivity models we obtained a universal high-f asymptote of the PSD. This behaviour is characterised by a 1/f² scaling, in analogy to Brownian motion [76] and scaled Brownian motion [78]. A first way to discriminate among models lies in the study of the ageing behaviour of the PSD, as already discussed in [78]. Indeed, we showed that the dependence of the PSD on the trajectory length T appears only for those random-diffusivity models that are characterised by an anomalous scaling of the MSD. We also showed that differences from one model to another appear in higher order moments. In particular, we obtained exact expressions for the coefficient of variation in all cases, proving that the latter can be a good indicator of the specific model (see table 1).

Finally, we established that the PDF of the random amplitude A carries most of the meaningful information. Namely, the coefficient of variation may be directly calculated from its moments. Moreover, its MGF is tightly related to that of the integrated diffusivity, thus reflecting the particular properties of the process Ψ_t. As we showed before [76,77], the distribution P(A) can even be evaluated meaningfully from experimental data of fairly short trajectories. In its useful role in data analysis the single-trajectory PSD approach thus complements other methods such as the time-averaged MSD and its amplitude variation [14,115,118–124], ageing analyses [77,115], or, in the context of non-Gaussian diffusion, the codifference methods [125].

The role of distinguishing different physical processes from measured single trajectory data therefore heavily lies on the amplitude fluctuations and the coefficient of variation encoded in them. This improved understanding of the single-trajectory PSD should therefore replace a common claim in many textbooks according to which the 1/f²-dependence of the spectrum was a fingerprint of Brownian motion. In line to previous works [76–78], in which we already alerted that this may be a deceptive concept, we have shown here that a wide range of random-diffusivity models with distinctly different behaviour and showing anomalous diffusion, exhibit precisely the 1/f²-dependence. Therefore, any experimental observation of the spectral density varying as 1/f² alone cannot be taken as proof of standard Brownian motion. One necessarily needs to consider the ageing behaviour of the spectral density, as in the case of superdiffusive fractional Brownian motion or scaled Brownian motion, evaluate the coefficient of variation of the spectral density, or determine the functional form of the PDF P(A) of the amplitude fluctuations.

In light of this the relevance of the presented results is twofold. First, they provide new and useful insights into the increasingly popular class of stochastic processes with random diffusivity used in the description of Fickian yet non-Gaussian diffusion in heterogeneous systems. Second, the results continue our ongoing analysis based on the single-trajectory PSD for different classes of stochastic processes, showing in particular the persistence of the 1/f²-scaling of the PSD, which appears to be robust—as long as we do not introduce correlations in the driving noise of the system, as studied in [77].

We note that if we redefine _t in the Langevin equation (3) as _t/S_t, we recover the seminal Black–Scholes (or Black–Scholes–Merton) equation for an asset price S_t with zero-constant trend and stochastic volatility $\sqrt{{D}_{0}{{\Psi}}_{t}}$ used in financial market models [126–128]. The relevance of diffusing-diffusivity approaches to economic and financial modelling was also discussed for the case of the squared Ornstein–Uhlenbeck process [61]. Namely, the resulting stochastic equation for D₀Ψ_t in this case is nothing else than the Heston model [83], a special class of the Cox–Ingersson–Ross model [82,129], and as such specifies the time evolution of the stochastic volatility of a given asset [63,83,130]. We note in this context that diffusing-diffusivity models are intimately related to random-coefficient autoregressive processes used in financial market analysis [131].

We finally note that realisation-to-realisation fluctuations of a stochastic process also turn out to be relevant in many scenarios of first-passage time statistics [132,133]. These fluctuations are connected to the typically broad ('defocused') PDFs of first-passage times even in simple geometries and the related feature of geometry-control [134–136]. It will be interesting to extend the existing first-passage time analyses of diffusing-diffusivity models [64,68] to the different processes studied herein, and to more complex geometries. In particular, given that for fractional Brownian motion the first-passage time distributions and the return-intervals distribution (the times between consecutive crossings of the diffusion process through a given level) have power-law exponents that are directly related to the scaling exponent of the power spectrum [137,138], one might speculate whether the corresponding distributions here are universal. However, the universal spectral behaviour predicted in this paper for a variety of random-diffusing models does not necessarily imply that the first-passage time properties will be the same for all the models. In fact, we expect that the first-passage time densities will exhibit different behaviours for the various models, to the same extent as the position PDFs unveiled here are significantly different from each other. The detailed analysis of the first-passage time properties of the new random-diffusivity models studied here constitutes one of the main directions of our future research.

VS is supported by the Basque Government through the BERC 2018-2021 program and by the Spanish Ministry of Economy and Competitiveness MINECO through BCAM Severo Ochoa accreditation SEV-2017-0718. FS acknowledges Padova University for support within grant PRD-BIRD191017. RM acknowledges Deutsche Forschungsgemeinschaft (DFG) for support within grant ME 1535/7-1, as well as the Foundation for Polish Science (Fundacja na rzecz Nauki Polskiej, FNP) for an Alexander von Humboldt Polish Honorary Research Scholarship. We acknowledge the support of the Deutsche Forschungsgemeinschaft (German Research Foundation) and the Open Access Publication Fund of Potsdam University.

It is convenient to rewrite formally the definition of the PSD in (8) in the form

$$\begin{equation}{S}_{T}\left(f\right)=\frac{1}{T}{\left({\int }_{0}^{T}\mathrm{d}t\;\mathrm{cos}\left(ft\right){x}_{t}\right)}^{2}+\frac{1}{T}{\left({\int }_{0}^{T}\mathrm{d}t\;\mathrm{sin}\left(ft\right){x}_{t}\right)}^{2}.\end{equation} \tag{ A.1 }$$

Our first step then consists in a standard linearisation of the expression in the exponential in (9). Taking advantage of the integral identity

$$\begin{equation}\mathrm{exp}\left(-{b}^{2}/\left(4c\right)\right)=\sqrt{\frac{c}{\pi }}{\int }_{-\infty }^{\infty }\mathrm{d}z\;\mathrm{exp}\left(-c{z}^{2}+\mathrm{i}bz\right)\end{equation} \tag{ A.2 }$$

for c > 0, we formally recast (9) into the form

$$\begin{equation}{\phi }_{\lambda }=\frac{1}{4\pi \lambda }{\int }_{-\infty }^{\infty }\mathrm{d}{z}_{1}{\int }_{-\infty }^{\infty }\mathrm{d}{z}_{2}\;\mathrm{exp}\left(-\frac{{z}_{1}^{2}+{z}_{2}^{2}}{4\lambda }\right){\left\langle \overline{\mathrm{exp}\left(\mathrm{i}{\int }_{0}^{T}\mathrm{d}t{Q}_{t}{x}_{t}\right)}\right\rangle }_{{\Psi}},\end{equation} \tag{ A.3 }$$

where Q_t is defined in (11). Now, the averaging over thermal noise realisations can be performed straightforwardly, yielding

$$\begin{align}\hfill \overline{\mathrm{exp}\left(\mathrm{i}{\int }_{0}^{T}\mathrm{d}t{Q}_{t}{x}_{t}\right)}& =\overline{\mathrm{exp}\left(\mathrm{i}\sqrt{2{D}_{0}}{\int }_{0}^{T}\mathrm{d}t\sqrt{{{\Psi}}_{t}}{\xi }_{t}{\int }_{t}^{T}\mathrm{d}\tau {Q}_{\tau }\right)}\hfill \\ \hfill & =\mathrm{exp}\left(-{D}_{0}{\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}{\left({\int }_{t}^{T}\mathrm{d}\tau {Q}_{\tau }\right)}^{2}\right),\hfill \end{align} \tag{ A.4 }$$

where we integrated by parts and used (3). Combining (A.3) and (A.4) we arrive at our general result (10).

The derivation of our main result (13) takes advantage of the explicit form of Q in (11) and the following calculation,

$$\begin{align}\hfill {\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}{\left({\int }_{t}^{T}\mathrm{d}\tau {Q}_{\tau }\right)}^{2}& =\frac{{z}_{1}^{2}}{{f}^{2}T}{\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}{\left(\mathrm{sin}\left(fT\right)-\mathrm{sin}\left(ft\right)\right)}^{2}+\frac{{z}_{2}^{2}}{{f}^{2}T}{\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}{\left(\mathrm{cos}\left(ft\right)-\mathrm{cos}\left(fT\right)\right)}^{2}\hfill \\ \hfill & \quad +\;\frac{2{z}_{1}{z}_{2}}{{f}^{2}T}{\int }_{0}^{T}\;\mathrm{d}t{{\Psi}}_{t}\left(\mathrm{sin}\left(fT\right)-\mathrm{sin}\left(ft\right)\right)\left(\mathrm{cos}\left(ft\right)-\mathrm{cos}\left(fT\right)\right).\hfill \end{align} \tag{ A.5 }$$

Inserting the latter expression into (10) and performing the integrations over z₁ and z₂ we find the expression in (13) with L_f(t₁, t₂) explicitly defined by

$$\begin{align}\hfill {L}_{f}\left({t}_{1},{t}_{2}\right)& =\frac{1}{2}\;\mathrm{cos}\left(2f{t}_{2}\right)-\frac{1}{2}\;\mathrm{cos}\left(2f{t}_{1}\right)-\frac{1}{2}\;\mathrm{cos}\left(f\left(T-{t}_{1}\right)\right)-\frac{1}{2}\;\mathrm{cos}\left(f\left(T-{t}_{2}\right)\right)-\frac{1}{4}\;\mathrm{cos}\left(2f\left(T-{t}_{1}\right)\right)\hfill \\ \hfill & \quad -\;\frac{1}{2}\;\mathrm{cos}\left(2f\left(T-{t}_{2}\right)\right)+\frac{1}{4}\;\mathrm{cos}\left(f\left(3T-{t}_{2}\right)\right)-\frac{1}{4}\;\mathrm{cos}\left(f\left(3T-{t}_{1}\right)\right)+\frac{3}{4}\;\mathrm{cos}\left(f\left(T+{t}_{1}\right)\right)-\frac{3}{4}\;\mathrm{cos}\left(f\left(T+{t}_{2}\right)\right)\hfill \\ \hfill & \quad -\;\frac{1}{2}\;\mathrm{cos}\left(f\left({t}_{1}-{t}_{2}\right)\right)-\frac{1}{4}\;\mathrm{cos}\left(2f\left({t}_{1}-{t}_{2}\right)\right)+\frac{1}{4}\;\mathrm{cos}\left(f\left(T-2{t}_{1}-{t}_{2}\right)\right)-\frac{1}{4}\;\mathrm{cos}\left(f\left(T-{t}_{1}-2{t}_{2}\right)\right)\hfill \\ \hfill & \quad +\frac{1}{2}\;\mathrm{cos}\left(f\left(T+{t}_{1}-2{t}_{2}\right)\right)+\frac{1}{2}\;\mathrm{cos}\left(f\left(T-2{t}_{1}+{t}_{2}\right)\right)+\frac{1}{4}\;\mathrm{cos}\left(f\left(T+2{t}_{1}-{t}_{2}\right)\right)\hfill \\ \hfill & \quad -\;\frac{1}{4}\;\mathrm{cos}\left(f\left(T-{t}_{1}+2{t}_{2}\right)\right)+\frac{1}{2}\;\mathrm{cos}\left(f\left(2T-{t}_{1}-{t}_{2}\right)\right)-\frac{1}{2}\;\mathrm{cos}\left(f\left(2T+{t}_{1}-{t}_{2}\right)\right)\hfill \\ \hfill & \quad +\;\frac{1}{2}\;\mathrm{cos}\left(f\left(2T-{t}_{1}+{t}_{2}\right)\right).\hfill \end{align} \tag{ A.6 }$$

For Riemann-integrable functions, according to the Riemann–Lebesgue lemma we have

$$\begin{equation}\underset{f\to \infty }{\mathrm{lim}}{\int }_{0}^{T}\mathrm{d}t{{\Psi}}_{t}\;\mathrm{cos}\left(f\left(T-t\right)\right)=0\end{equation} \tag{ B.1 }$$

with probability 1 (however, nothing can be said about how fast zero is approached in the general case). Once (B.1) holds, one finds that for L_f(t₁, t₂) defined in (A.6),

$$\begin{equation}\underset{f\to \infty }{\mathrm{lim}}{\int }_{0}^{T}\mathrm{d}{t}_{1}{{\Psi}}_{{t}_{1}}{\int }_{0}^{T}\mathrm{d}{t}_{2}{{\Psi}}_{{t}_{2}}{L}_{f}\left({t}_{1},{t}_{2}\right)=0.\end{equation} \tag{ B.2 }$$

Our expression for the PDF P(A) in (26) and (27) rely on the following series expansion of the product of two Bessel functions,

$$\begin{equation}{J}_{0}\left(\left(1+\frac{1}{\sqrt{3}}\right)\sqrt{2zA}\right){J}_{0}\left(\left(1-\frac{1}{\sqrt{3}}\right)\sqrt{2zA}\right)=\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^{n}}{{\left(n!\right)}^{2}}{\left(\frac{\sqrt{3}+1}{\sqrt{6}}\right)}^{2n}\;{}_{2}F_{1}\left(-n,-n;1;\frac{1-\sqrt{3}/2}{1+\sqrt{3}/2}\right){\left(zA\right)}^{n}.\end{equation} \tag{ C.1 }$$

The form of the PDF in (95) stems from the expansion

$$\begin{equation}\frac{1}{\sqrt{\mathrm{cosh}\left(\sqrt{4{D}_{\mathrm{B}}Tz/{a}^{2}}\right)}}=\sqrt{2}\sum _{n=0}^{\infty }\left(\genfrac{}{}{0pt}{}{-1/2}{n}\right)\;\mathrm{exp}\left(-2\sqrt{\frac{{D}_{\mathrm{B}}Tz}{{a}^{2}}}\left(2n+\frac{1}{2}\right)\right).\end{equation} \tag{ C.2 }$$

The result in (69) involves toroidal functions defined by

$$\begin{equation}{P}_{n-1/2}\left(\mathrm{cosh}\left(\eta \right)\right)=\frac{2}{\pi }{\mathrm{e}}^{-\left(n+1/2\right)\eta }{\int }_{0}^{\pi /2}\frac{\mathrm{d}\phi }{{\left(1-2{\mathrm{e}}^{-\eta }\;\mathrm{sinh}\left(\eta \right)\;{\mathrm{sin}}^{2}\left(\phi \right)\right)}^{n+1/2}}.\end{equation} \tag{ C.3 }$$

Setting exp(−η) sinh(η) = 1/3, i.e., $\eta =\mathrm{ln}\left(\sqrt{3}\right)$, we obtain expression (69).