Universal spectral features of different classes of random diffusivity processes

Stochastic models based on random diffusivities, such as the diffusing-diffusivity approach, are popular concepts for the description of non-Gaussian diffusion in heterogeneous media. Studies of these models typically focus on the moments and the displacement probability density function. Here we develop the complementary power spectral description for a broad class of random diffusivity processes. In our approach we cater for typical single particle tracking data in which a small number of trajectories with finite duration are garnered. Apart from the diffusing-diffusivity model we study a range of previously unconsidered random diffusivity processes, for which we obtain exact forms of the probability density function. These new processes are different versions of jump processes as well as functionals of Brownian motion. The resulting behavior subtly depends on the specific model details. Thus, the central part of the probability density function may be Gaussian or non-Gaussian, and the tails may assume Gaussian, exponential, log-normal or even power-law forms. For all these models we derive analytically the moment-generating function for the single-trajectory power spectral density. We establish the generic $1/f^2$-scaling of the power spectral density as function of frequency in all cases. Moreover, we establish the probability density for the amplitudes of the random power spectral density of individual trajectories. The latter functions reflect the very specific properties of the different random diffusivity models considered here. Our exact results are in excellent agreement with extensive numerical simulations.

Stochastic models based on random diffusivities, such as the diffusing-diffusivity approach, are popular concepts for the description of non-Gaussian diffusion in heterogeneous media. Studies of these models typically focus on the moments and the displacement probability density function. Here we develop the complementary power spectral description for a broad class of random diffusivity processes. In our approach we cater for typical single particle tracking data in which a small number of trajectories with finite duration are garnered. Apart from the diffusing-diffusivity model we study a range of previously unconsidered random diffusivity processes, for which we obtain exact forms of the probability density function. These new processes are different versions of jump processes as well as functionals of Brownian motion. The resulting behavior subtly depends on the specific model details. Thus, the central part of the probability density function may be Gaussian or non-Gaussian, and the tails may assume Gaussian, exponential, log-normal or even power-law forms. For all these models we derive analytically the moment-generating function for the single-trajectory power spectral density. We establish the generic 1/f 2 -scaling of the power spectral density as function of frequency in all cases. Moreover, we establish the probability density for the amplitudes of the random power spectral density of individual trajectories. The latter functions reflect the very specific properties of the different random diffusivity models considered here. Our exact results are in excellent agreement with extensive numerical simulations.

I. INTRODUCTION
Diffusive processes came to the attention of the broader scientific community with the experiments on "active molecules" by Robert 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 Adolf 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 Karl 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 molecular tracers in complex environments such as living biological cells [8][9][10][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 phe-nomena in the stochastic particle motion and thus call for new theoretical concepts [14][15][16] on top of already known approaches [17][18][19].
Among these new insights is that endogenous and introduced tracers in living biological cells perform anomalous diffusion of the form r 2 (t) K α t α 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][23][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 nanometer 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], or 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" 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][42][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 behaviors from a non-Gaussian shape at shorter time scales to an effective Gaussian behavior 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][50][51][52]54]. A microscopic realization of such a behavior was proposed for a model of diffusion during a polymerization process [53]. However, in superstatistics (and in the related process called generalized grey Brownian motion [55][56][57]) the distribution is a constant of the motion and thus no crossover behavior 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 characterized 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 analyzed further, and their application is mainly the diffusive dynamics in heterogeneous systems [59][60][61][62][63][64][65][66][67][68][69][70].
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 precide 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 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 [71][72][73][74][75][76], 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 (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. 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) as functions of frequency f and measurement time T . We previously analyzed the behavior of S(f, T ) for different diffusion scenarios [71][72][73] and demonstrated that it is practically useful in the analysis of experimental data [71,72]. In what follows we derive the momentgenerating 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 II with a description of the model and in section III we report general results on the spectral properties of this class of processes. Specific examples of diffusingdiffusivity models are described in sections IV-VI. The first example is the well known case in which the diffusivity is modeled as the squared Ornstein-Uhlenbeck process. In the second group of examples we analyze 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 VII we draw our conclusions. In the appendix we report details on the explicit derivations of our results.

II. RANDOM DIFFUSIVITY PROCESSES
We consider a class of one-dimensional stochastic processes x t that obey the Langevin equation in the Itô convention,ẋ Here D 0 is a constant, dimensional coefficient in units length 2 /time, and in our analysis we will assume the initial condition x 0 = 0. In equation (3) ξ t denotes a standard Gaussian white noise with zero mean and covariance function ξ t ξ t = δ (t − t ). The bar here and henceforth denotes averaging with respect to the noise ξ t . Lastly Ψ t is a positive-definite random function, which multiplies D 0 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 T 0 dtΨ 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 α produces so-called scaled Brownian motion [77,78]. 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 jumpprocess, that attains independent, identically distributed random values. We also present several examples when Ψ t is subordinated to standard unbiased Brownian mo- 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 realization of the noise and a given realization of Ψ t , to obtain The characteristic function of x t can be written down in the form where the bar stands for averaging over thermal histories, while the angular brackets denote averaging over the realizations of the random function Ψ t . The thermal average can be performed straightforwardly to give The desired PDF Π(x, t) can then be written as In the following section IV 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.

III. GENERAL THEORY
We first obtain exact expressions for the PSD (2) and then study the limiting behavior for high frequencies.
A. 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 ), 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 with λ ≥ 0. Relegating some intermediate calculations to Appendix A we find the following expression for φ λ in (9) averaged over thermal noises, where Performing the inverse Laplace transform of expression (10), we find the general result for the PDF where J 0 (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) where L f (t 1 , t 2 ) is defined by the somewhat lengthy expression (A6). 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 realization 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: in particular, for Ψ t = 1 it describes the MGF in case of standard Brownian motion [71], while the choice Ψ t = t α−1 corresponds to the case of scaled Brownian motion recently studied in [73].

B. Exact high frequency limiting behavior
As already remarked we here concentrate on random processes Ψ t which, for any finite T , are Riemannintegrable 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 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 0 λ/f 2 so that the large-f spectrum of a single-trajectory PSD has the universal form regardless of the specific choice of Ψ t . Here A is a dimensionless, random amplitude, which differs from realization to realization. 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 [71][72][73].
The MGF Φ λ of the random amplitude A follows from (14) and can be written as where I 0 (z) is the modified Bessel function of the first kind, and is the MGF of the integrated diffusivity (see [64,66]) Relation (16) tightly links the MGFs of A and τ T . Moreover, note that the MGF of the diffusing-diffusivity model in (5) is tightly related to the one in (17): specifically, 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 highfrequency behavior of the PSD. Taking the inverse Laplace transform with respect to the parameter λ we evaluate the PDF of A, This exact expression determines the large-f behavior of the PDF p(S T (f ) = S), as f → ∞. The exact large-f forms in (16) and (20) will serve as the basis of our analysis for several particular choices of the process Ψ t in section IV. Before proceeding, we stop to make several general statements.
(i) 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, where 2 F 1 (a, b; c; z) is the Gauss hypergeometric function.
Since the moments of τ T are related to the moments of the process x T [66] we also find (ii) Equations (15) and (22) permit us to directly access the coefficient of variation γ of the PDF p(S T (f ) = S) in the large-f limit. We get straightforwardly 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 0 and f when f is only large enough.
(iii) The behavior 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 (C1) in Appendix C). Inserting the expansion in (C1) in (20) and integrating over z we find 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 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 (26).

IV. DIFFUSIVITY MODELED AS SQUARED ORNSTEIN-UHLENBECK PROCESS
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 highfrequency 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. We note that diffusing-diffusivity models are intimately related to random-coefficient autoregressive processes [79]. The 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 Ψ t = Y 2 t is one of the most common models of diffusing diffusivity, which satisfies, due to the Itô's formula, where τ = τ /2, σ = √ 2σ , andΨ = σ 2 τ /2 = σ 2 τ /2. This model was extended in [59][60][61] by considering Ψ t as the sum of n independent squared Ornstein-Uhlenbeck processes, when (28) still holds withΨ = nσ 2 τ /2. More generally, settingΨ to be any positive constant, the Langevin equation (28) defines the so-called Feller process [80], also known as square root process or the Cox-Ingersoll-Ross process [81]. 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 [82] and employed in [59][60][61]. Its computation for the Feller process in (28) was presented in [63], where ω = √ 1 + 4σ 2 τ 2 λ and ν =Ψ/(τ σ 2 ). In particular, settingΨ = σ 2 τ /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 x 2 T and x 4 T from [61,63] we get from (24) that 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 Fig. 1, in which we observe excellent agreement between the simulations and the theoretical results. The 1/f 2 scaling is recovered for any value of τ . The coefficient of variation γ converges to different values when we change τ , according to (30). Note that this result reflects the different degrees of broadness of the PDF of the random amplitude A.

V. DIFFUSIVITY MODELED AS A JUMP PROCESS
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 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 realization of the process Ψ t we thus have and hence Evaluating explicitly the derivatives ∂ λ Υ(T ; λ) and ∂ 2 λ Υ(T ; λ) at λ = 0, we get when the first two moments E{ψ k } and E{ψ 2 k } exist. Let us consider two representative cases.

Example I: Gamma distribution
First, we consider the Gamma distribution, with the shape parameter ν > 0 and the scale parameter ψ 0 > 0. From (33), we deduce 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 0 w 2 ), we get where the normalization coefficient is With the properties for |z| → 0 and ν > 1, as well as for |z| → ∞, the asymptotic behaviors of the PDF are given by for |x| → 0 and νt > 3δ/2, as well as for |x| → ∞.
The functional behavior of the PDF Π(x, t) is shown in Fig. 2. We see that by changing δ we can observe different shapes of Π(x, t). When δ = 1 (panels a) and c)) the Gaussian approximation (49) 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 behaviors at small and large x are well described by the asymptotic expansions in (41) and (42).

Amplitude-PDF P (A)
The MGF of the amplitude A of the jump processdiffusivity model is given by so that In particular one has E{ψ k } = νΨ 0 and E{ψ 2 k } = Ψ 2 0 ν(ν + 1) and thus In the limit δ → 0 and N → ∞, with δN = T fixed, we have Hence, and  (41) and (42).
This means that we have essentially the same behavior as for standard one-dimensional Brownian motion, however, with renormalized coefficients (compare with the result in [71]). Indeed, if we use (46) and recall that Φ w = Υ(T ; D 0 w 2 ), we readily obtain In Fig. 3 we show a direct comparison between the numerical and theoretical results for the Gamma distribution with ψ 0 = 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 (35), while the mean behavior is not.

Example II: Lévy-Smirnov distribution
In our second example we consider the Lévy-Smirnov distribution for which Eq. (33) yields Position-PDF Π(x, t) As a consequence, we obtain the following analytical expression for the PDF, where we recognize the power-law behavior, that is already built into relation (50). Note that expression (52) represents the Cauchy distribution, whose median grows with time t. The PDF is shown for two different values of δ in Fig. 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 behavior (52) is readily discernible.

Amplitude-PDF P (A)
The MGF for the random amplitude A reads and with ξ = (4Aδ)/(3ψ 0 T ). Note that in the limit A → ∞, the leading behavior of P (A) follows Thus, the PDF P (A) inherits the property of diverging moments from the parental Lévy-Smirnov distribution. Figure 5 summarizes the properties of the PSD for the jump process with Lévy-Smirnov distribution (ψ 0 = 0.5). We observe that, despite the fat-tailed PDF in (52) we still observe the universal 1/f 2 scaling of the PDF. Concurrently, the PDF of the random amplitude A features the power-law behavior according to (54). Note that the non-existence of the moments of P (A) generates a pronounced scatter in the amplitude of the average power spectrum.

VI. DIFFUSIVITY MODELED AS A FUNCTIONAL OF BROWNIAN MOTION
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 We where V is some prescribed, positive-definite function. Note that random variables of the form T 0 dtV [B t ] appear across many disciplines, including probability theory, statistical analysis, computer science, mathematical finance and physics. Starting from earlier works [83,84,86,87,89], 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., [90][91][92][93][94][95] 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.
Example I: Ψt = B 2 t /a 2 First, we consider squared Brownian motion where a is a parameter of unit length. Note that here the process x t in (3) is superdffusive, such that its meansquared displacement obeys 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 [84,85] (see also [86]). In our notation, and the PDF Π(x, t) for this process is given by where c = 2 √ D 0 D B /a and P ν (z) is the Legendre function of the first kind. The latter admits the representa- In the limits, for |x| → 0, and for |x| → ∞. As a consequence, the behavior of the PDF for small x is approximately Gaussian, exp(−const.x 2 ), where const. is a constant that can be expressed via the polylogarithm function. Conversely, for |x| → ∞. Note that the integrals exist for any n > 0 and, hence, all negative moments of T 0 dtB 2 t exist, as well. As a consequence, P (A) is an analytic function of A. We immediately obtain the coefficient of variation from (24) as γ = √ 17/2. The shape of the PDF is shown in Fig. 6. We clearly observe that the central part is approximately Gaussian while the tails have an exponential trend, as expressed explicitly by the asymptote (65).

Amplitude-PDF P (A)
Using (59) we then find that the MGF of the random amplitude A and the corresponding PDF are given by the integrals and The series representation of P (A) in (68) can be found directly by taking advantage of expansion (C2) in Appendix C and our result (54), to yield with ξ n = a 2 A Example II: Ψt = Θ(Bt) As the second example we choose the cut-off Brownian motion 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 T 0 dtΨ 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. The MSD of the process x t , as one can straightforwardly check, is just that is, a standard diffusion law in which the diffusion coefficient is reduced by the factor 1/2. For higher order moments one expects, of course, significant departures from the standard diffusive behavior.
The MGF of the random variable τ T = T 0 dtΘ(B t ), which has a bounded support on (0, T ), has been first derived by Kac [86], and Erdös and Kac [87]. Rewriting their result in our notation we have Note that the inverse Laplace transform of this expression produces the celebrated Lévy arcsine law [88]. With result (73), the desired PDF is given by Recalling that K 0 (z) ∼ − ln(z/2) − γ EM for |z| → 0, where γ EM is the Euler-Mascheroni constant, for small x we have For large x we use (40) and obtain the asymptotic behavior for the PDF, Note that the integrals ∞ 0 dλλ n Υ(T ; λ) diverge for any n > 0, which means that τ T does not have negative moments. One therefore expects that P (A) is a nonanalytic, diverging function in the limit A → 0.
The PDF for this process is shown in Fig. 8. We see that the central part of the PDF is strongly non-Gaussian, while the tails are Gaussian, in agreement with the asymptotic behaviors (75) and (76).

Amplitude-PDF P (A)
Inserting expression (73) into (16) and performing the integration over z, we arrive at the following, remarkably compact expression for the MGF, where K(x) is the complete elliptic integral of the first kind, Note that the large-f asymptotic form in (77) 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 (77). After some formal manipulations this yields the following expression for the PDF, Multiplying both sides of the latter equation by A n and integrating over A from 0 to ∞, we get the following simple expression for the moments of the random amplitude A of an arbitrary order, Then, from (24), we readily get the coefficient of variation, γ = 19/8. The small-A asymptotic behavior of P (A) can be deduced directly from (79). Expanding the exponential function in the integral into the Taylor series in powers of A and expressing the emerging generalized elliptic integrals via their representations in terms of the toroidal functions P n−1/2 (cosh(η)) (see (C3) in Appendix C), we get (81) For the opposite limit A → ∞, we conveniently rewrite Eq. (54) in the form where L (−1/2) n (x) are associated Laguerre polynomials. We focus next on the asymptotic behavior of the function g(u) = u −1/2 ∞ n=0 Γ(n + 1/2) n! = 1 s 1/2 ∞ n=0 Γ 2 (n + 1/2) (n!) 2 2(s − 1) 3s where K is the complete elliptic integral defined in Eq. (78). In the limit s → 0 (corresponding to A → ∞), .
(85) Inverting the Laplace transform and integrating over φ, we find Thus, in the limit A → ∞, the leading behavior of the PDF P (A) yields in the form Figure 9 summarizes the numerical results for this case. Again we observe excellent agreement with the theoretical results.
Finally we link the process Ψ t to so-called geometric Brownian motion and set where a is a parameter of unit length. In this case, x t exhibits an anomalously strong superdiffusion such that Position-PDF Π(x, t) Random variables of this form have been widely studied in the mathematical finance literature (see, e.g., [90]). 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 [96][97][98] (see also [99,100]) and also appears in different contexts in the analysis of transport phenomena in disordered media (see, e.g., [101][102][103][104][105][106][107]) as well as characterizes some features of the melting transition of heteropolymers [108]. The Laplace transform of the time-averaged geometric Brownian motion in our notation reads [102][103][104] where K ix is the modified Bessel function of the second kind with purely imaginary index. As a consequence, the PDF is given by where Recalling that arcsinh(z) ∼ z for z → 0 and arcsinh(z) ∼ ln(2z) for z → ∞, we express the asymptotic behavior of the PDF as for |x| → 0 and for |x| → ∞.
The PDF is shown in Fig. 10. In this case, according to the asymptotic expansions (93) and (94) we observe that the central part of the PDF is approximately Gaussian, while the tails follow a log-normal shape.

VII. CONCLUSIONS
Quite typically stochastic time series are evaluated in terms of the ensemble averaged MSD x 2 t . 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 realized that the amplitude fluctuations between individual trajectories in fact harbors important quantitative information characteristic for a given stochastic process [14,15,[109][110][111].
Similar to the consideration of time averaged MSDs for trajectories of finite measurement time T , we here analyzed the single-trajectory PSD of stochastic trajectories characterized by random diffusivities. Following our previous work on standard Brownian motion [71], as well as fractional [72] and scaled Brownian motion [73], we here investigated the detailed behavior of singletrajectory 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 describe a general procedure to obtain the PDF Π(x, t) for all such models, showing that different choices of the underlying diffusivity Ψ t lead to very different resulting behaviors.
We started our discussion from the by-now wellestablished and widely studied model of "diffusingdiffusivity", 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 modeled 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 large-f asymptote of the PSD. This behavior is characterized by a 1/f 2 scaling, in analogy to Brownian motion [71] and scaled Brownian motion [73]. 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.
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 [71,72], the distribution P (A) can even be evaluated meaningfully from experimental data of fairly short trajectories. In its useful role in data analysis the singletrajectory PSD approach thus complements other methods such as the time-averaged MSD and its amplitude variation [14,109,[112][113][114][115][116][117][118], ageing analyses [72,109], or, in the context of non-Gaussian diffusion, the codifference methods [119].
The role of distinguishing different physical processes from measured single trajectory data therefore heaviliy lies on the amplitiude 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 2 -dependence of the spectrum was a fingerprint of Brownian motion. In line to previous works [71][72][73], 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 behavior and showing anomalous diffusion, exhibit precisely the 1/f 2 -dependence. Therefore, any experimental observation of the spectral density varying as 1/f 2 alone cannot be taken as proof of standard Brownian motion. One necessarily needs to consider the ageing behavior 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 2 -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 [72].
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 √ D 0 Ψ t used in financial market models [120][121][122]. The relevance of diffusing-diffusivity approaches to economic and fincancial modeling was also discussed for the case of the squared Ornstein-Uhlenbeck process [61]. Namely, the resulting stochastic equation for D 0 Ψ t in this case is nothing else than the Heston model [123], a special class of the Cox-Ingersson-Ross model [124,125], and as such specifies the time evolution of the stochastic volatility of a given asset [63,123,126].
We finally note that realization-to-realization fluctuations of a stochastic process also turn out to be relevant in many scenarios of first-passage time statistics [127,128]. These fluctuations are connected to the typically broad ("defocused") PDFs of first-passage times even in simple geometries [129,130]. 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. (A1) Our first step then consists in a standard linearization of the expression in the exponential in (9). Taking advantage of the integral identity for c > 0, we formally recast (9) into the form where Q t is defined in (11). Now, the averaging over thermal noise realizations can be performed straightforwardly, yielding where we integrated by parts and used (3). Combining (A4) and (A3) we arrive at our general result (10).