Power spectral density of a single Brownian trajectory: what one can and cannot learn from it

The power spectral density (PSD) of a stochastic process X_t, which is formally defined as (see, e.g., [1])

$$\begin{eqnarray}&&{\mu }_{S}(f,\infty )=\mathop{\mathrm{lim}}\limits_{T\to \infty }\displaystyle \frac{1}{T}{\mathbb{E}}\left\{{\left|{\int }_{0}^{T}{{\rm{e}}}^{{\rm{i}}{f}{t}}{X}_{t}{\rm{d}}{t}\right|}^{2}\right\},\end{eqnarray} \tag{ 1 }$$

provides important insights into the spectral content of X_t. Here and in what follows, the symbol ${\mathbb{E}}\{\ldots \}$ denotes averaging with respect to all possible realizations of the process, i.e., the expectation. Often the PSD as defined in equation (1) has the form ${\mu }_{S}(f,\infty )=A/{f}^{\beta }$, where A is an amplitude and β, (typically, one has 0 < β ≤ 2 [2]), is the exponent characteristic of the statistical properties of X_t. In experiments and numerical modeling, the PSD is determined using a periodogram estimate for a wide variety of systems in physics, biophysics, geology etc. A by far non-exhaustive list of systems exhibiting 1/f^β spectra includes electrical signals in vacuum tubes, semiconductor devices, metal films and other condensed matter systems [3–5], quantum dots [6–8], nano-electrodes [9] and two-dimensional graphene layers with a widely tunable concentration of carriers [10]. The PSD has been analyzed, as well, for the trajectories of tracers in artificial crowded fluids [11], for active micro-rheology of colloidal suspensions [12], Kardar–Parisi–Zhang interface fluctuations [13], for sequences of earthquakes [14], weather data [15–17], biological evolution [18], human cognition [19], network traffic [20] and even for the loudness of music recordings (see, e.g., [21, 22]).

On the theoretical side, the PSDs showing the 1/f^β dependences have been calculated analytically for diverse situations, including, e.g., the dynamics in chaotic Hamiltonian systems [23], periodically-driven bistable systems [24], fluctuations of a phase-separating interface [25], several diffusive and non-diffusive transport processes (see, e.g., [8, 26–29]), the running maximum of Brownian motion (BM) [30], diffusion in presence of strong quenched disorder [31–34] and the electric-field-driven ion transport through nanometer-scale membrane pores [35]. The PSD of BM in an optical trap has been scrutinized in [36, 37], which permitted the calibration of optical tweezers, making them a powerful tool for force spectroscopy, local viscometry, and other applications (see, e.g. [38, 39]).

The systematic measurements of single particle trajectories started with Perrin more than a century ago [40]. Just a few years later Nordlund [41] developed impressive experimental techniques, followed by a line of further refinements up to Kappler's measurements [42], to record sufficiently long trajectories to enable quantitative analysis on the basis of individual trajectories—without the need of prior ensemble averaging. Nowadays, with the advent of modern microscopy and supercomputing, scientists routinely measure long trajectories X_t of submicron tracer particles or even single molecules [11, 39, 43–47].

In parallel to this experimental progress, there was a general shift of interest towards understanding statistical properties of individual realizations of stochastic processes. In particular, a conceptually important question often raised within this context is if one can reliably extract information about the ensemble-averaged properties of random processes from single-trajectory data [46]. Considerable theoretical progress has been achieved, for instance, in finding the way to get the ensemble-averaged diffusion coefficient from a single Brownian trajectory, which task amounts to seeking properly defined functionals of the trajectory which possess an ergodic property (see, e.g., [43–46, 48–56] for a general overview).

One of such functionals used in single-trajectory analysis is the time-averaged mean squared displacement (MSD) (see, e.g., [39, 43–47])

$$\begin{eqnarray}&&\overline{{\delta }^{2}(\tau )}=\displaystyle \frac{1}{T-\tau }{\int }_{0}^{T-\tau }{({X}_{t+\tau }-{X}_{t})}^{2}{\rm{d}}{t},\end{eqnarray} \tag{ 2 }$$

where the bar denotes time averaging and the lag time τ acts as the size of the analysis window sliding over the time series X_t. For finite observation time T, the time-averaged MSD varies randomly for different trajectories even under identical physical conditions. For both normal and anomalous diffusion this variation is mainly seen as an amplitude scatter of $\overline{{\delta }^{2}(\tau )}$ at given lag times, which remains remarkably constant (apart of relatively small local fluctuations) as function of τ [57–60].

For the time-averaged MSD the fluctuations are often quantified in terms of the so-called ergodicity breaking parameter $\mathrm{EB}={\mathbb{E}}\{{\xi }^{2}\}-1$, where the dimensionless variable $\xi =\overline{{\delta }^{2}(\tau )}/{\mathbb{E}}\{\overline{{\delta }^{2}(\tau )}\}$ for a given τ measures the deviations from the trajectory-average ${\mathbb{E}}\{\overline{{\delta }^{2}(\tau )}\}=(1/N){\sum }_{i=1}^{N}\overline{{\delta }_{i}^{2}(\tau )}$ [57, 61], with N being the number of observed trajectories. Clearly, $\mathrm{EB}\equiv {\gamma }^{2}$ with γ being the coefficient of variation of the distribution ${P}(\bar{{{\delta }}^{2}({\tau })})$.

For BM, the fluctuations decrease with growing observation time such that $\mathrm{EB}\sim 4\tau /(3T)$ [44, 50, 57], which permits to deduce the diffusion coefficient specific to an ensemble of trajectories from just a single, sufficiently long, trajectory. A similar decay of EB to zero for large T is observed for fractional BM [62] and scaled BM [63]. However, for anomalous diffusion with scale-free distributions of waiting times [61] and processes with systematic, spatially varying diffusion coefficient [64] the inequality $\mathrm{EB}\ne 0$ persists even in the limit $T\to \infty$, which reveals ageing in the sense that the properties of the system, such as the effective diffusivity, are perpetually changing during the measurement and depend on T [46, 57, 65]. In these systems, the fluctuations as measured by EB do not decay to zero with increasing T, and thus time averages of physical observables remain random quantities, albeit with well defined distributions (see, e.g. [50, 57, 61]). Many other facets of the time-averaged MSD in equation (2) and of the parameter EB were exhaustively studied to date for a variety of diffusive processes, providing a solid mathematical framework for the analysis of individual trajectories [39, 57, 62, 66–69].

Single-trajectory PSDs have been studied in several cases. In particular, the PSD of the loudness of musical recordings was discussed in [21, 22] and the spectra of temperature data were presented in [16, 17]. We note that, of course, for these two examples averaging over 'an ensemble' does not make any sense since the latter simply does not exist, likewise, e.g., the case of the financial markets data for which the single-trajectory MSD has been recently analyzed [70]. Further on, power spectra of individual time-series were examined for a two-state stochastic model describing blinking quantum dots [26] and also for single-particle tracking experiments with tracers in artificially crowded fluids [11]. It was shown that, surprisingly enough, the estimation of the exponent characterizing the power spectrum of an ensemble from the single-trajectory PSD is rather robust. Next, in [71, 72] the power spectra of the velocities of independent motile amoeba (see also [47, 73]) were analyzed, revealing a robust high-frequency asymptotic 1/f² form, persisting for all recorded individual trajectories.

Clearly, these numerical observations raise several challenging questions. Indeed, could it be possible that the exponent characterizing the frequency-dependence of the standard PSD, defined as an ensemble-average property, can be observed already on a single-trajectory level? If true, does it hold for any stochastic transport process or just for some particular examples? Moreover, in which range of frequencies can such a behavior be observed? What are the distributions of the amplitudes entering the relation between a single-trajectory PSD and its ensemble-average counterpart and how broad are they? Evidently, numerical analyzes may give a certain degree of understanding of some particular features, but a deep insight can be obtained only in conjunction with a full mathematical solution. Unlike the single-trajectory MSD, for which a deep knowledge has been already acquired via an exact analytical analysis, a similar analysis of the statistical properties of a single-trajectory PSD is lacking at present, although its ensemble-averaged counterpart is widely used as an important quantifier of different properties of random trajectories in diverse areas of engineering, physics and chemistry.

In this paper, going beyond the text-book definition (1), we concentrate on the question what information can be reliably obtained if one defines the PSD of a single, finite-time realization of X_t. We here focus on the paradigmatic process of BM. This choice is two-fold: first, BM is ubiquitous in nature which renders this analysis particularly important. Second, it permits us to obtain an exact mathematical solution of the problem: we calculate exactly the moment-generating function and the full probability density function of the single-trajectory PSD and its moments of arbitrary order in the most general case of arbitrary frequency, arbitrary (finite) observation time and arbitrary number of the projections of a d-dimensional BM onto the coordinate axes. This furnishes yet another example of a time-averaged functional of BM, whose moment-generating function and full probability distribution can be calculated exactly (see, e.g., [74]).

Capitalizing on these results, we observe that for a sufficiently large T (and frequency f bounded away from zero) for any realization of the process a single-trajectory PSD is proportional to its first moment (1), and the latter embodies the full dependence on the frequency and on the diffusion coefficient specific to an ensemble of the trajectories. This means that the frequency-dependence can be deduced from a single trajectory. In other words, there is no need to perform averaging over an ensemble of trajectories—one long trajectory suffices. However, the proportionality factor, connecting a single-trajectory PSD and its ensemble-averaged counterpart—a numerical amplitude—is random and varies from realization to realization. Due to this fact, one cannot infer the value of the ensemble-averaged diffusion coefficient from the amplitude of a single-trajectory PSD. The distribution function of this amplitude is calculated exactly here and its effective width is quantified using standard criteria.

As a proof of concept, we revisit our predictions for a continuous-time BM—an idealized process with infinitesimal increments—resorting to a numerical analysis based on Monte Carlo simulations of discrete-time random walks, and also using experimental single-trajectory data for the diffusive motion of micron-sized polystyrene beads in a flow cell. We demonstrate that our theoretical prediction for the relation connecting a single-trajectory PSD and its ensemble-averaged counterpart is corroborated by numerical and experimental results. Additionally we show that the predicted distribution of the amplitude is consistent with numerical results, which means that the framework developed here is justified and allows for a meaningful analysis of experimental trajectories.

Pursuing this issue further, we address several general questions emerging in connection with a comparison of our analytical predictions against numerical simulations and experimental data. To this end, we first consider Wiener's representation of BM in form of an infinite Fourier series with random coefficients, whose truncated version is often used in numerical simulations. We show that the distribution of the single-trajectory PSD obtained from Wiener's representation in which just N terms are kept, instead of an infinite number, has exactly the same form as the one obtained for the continuous-time BM when fT is within the interval (0, πN). Outside of this interval, the probability density function of the truncated PSD converges to a different form.

We examine the case when a trajectory of a continuous-time BM is recorded at some discrete time moments, so that it is represented by a set of M points. The single-trajectory PSD, i.e., the periodogram, becomes a periodic function of fT with the prime period equal to 2πM. We analyze several aspects of this discrete-time problem: we study how large M should be taken at a fixed observation time T so that we may recover the results obtained for the continuous-time BM. We analyze the limiting forms of the distribution of a single-trajectory periodogram and show, in particular, that for fT kept fixed and $M\to \infty$, the distribution converges to the form obtained for a continuous-time BM. On the contrary, when fT is left arbitrary so that it may assume any value within the prime period, that is, ${fT}\in (0,2\pi M)$, the distribution of a single-trajectory periodogram converges to a different limiting form as $M\to \infty$. Our analysis demonstrates that when fT belongs to a certain interval within the prime period, a single-trajectory periodogram equals, up to a random numerical amplitude, the ensemble-averaged periodogram, and the latter embodies the full dependence on f and T. Therefore, similarly to the continuous-time case, the correct spectrum can be obtained already from a single trajectory. These new results are the foundation for the use of single trajectory PSD in the quantitative analysis of single or few recorded particle trajectories, complementing the widely used concept of the single trajectory MSD.

The paper is outlined as follows: in section 2 we introduce our basic notations. In section 3 we first derive explicit expressions for the variance of a single-trajectory PSD and the corresponding coefficient of variation of the probability density function, and present exact results for the moment-generating function of a single-trajectory PSD, its full probability density function and moments of arbitrary order (section 3.1). Section 3.2 is devoted to the relation connecting a single-trajectory PSD and its ensemble-averaged counterpart, while section 3.3 discusses fluctuations of the amplitude in this relation. Next, in section 4 we analyze the probability density function of a single-trajectory PSD obtained by truncating Wiener's representation of a continuous-time BM. Section 5 presents an analogous analysis for the case when a continuous-time trajectory is recorded at discrete time moments. In section 6 we conclude with a brief recapitulation of our results and outline some perspectives for further research. Additional details are relegated to appendix A, in which we present exact results for the distributions in the special case f = 0, and to appendix B, where we discuss several cases in which the full distribution of a single-trajectory periodogram can be evaluated exactly in the discrete-time settings.

Let X_t, t ∈ (0, T) denote a Brownian trajectory in a d-dimensional continuum and ${X}_{t}^{(j)}$ with j = 1, 2, ..., d stand for the projection of ${{\bf{X}}}_{t}$ on the axis x_j. The projections ${X}_{t}^{(j)}$ are statistically independent of each other and (likewise the BM X_t itself) each projection ${X}_{t}^{(j)}$ is a Gaussian process with zero mean and variance

$$\begin{eqnarray}&&{\mathbb{E}}\{{({X}_{t}^{(j)})}^{2}\}=2{Dt},\end{eqnarray} \tag{ 3 }$$

where D is the diffusion coefficient of ${{\bf{X}}}_{t}$. In a random walk sense, $D=\langle \delta {x}^{2}\rangle /(2d\langle \delta \tau \rangle )$, where $\langle \delta {x}^{2}\rangle$ is the variance of the step length and $\langle \delta \tau \rangle$ is the mean time interval between consecutive steps. Hence, D contains the dimension d of the unprojected motion.

In text-book notations, the PSD of each of the projections is defined by equation (1), that is,

$$\begin{eqnarray}&&{\mu }_{S}^{(j)}(f,\infty )=\mathop{\mathrm{lim}}\limits_{T\to \infty }{\mu }_{S}^{(j)}(f,T),\end{eqnarray} \tag{ 4 }$$

where (taking into account that ${X}_{t}^{(j)}$ is real-valued) the T-dependent function ${\mu }_{S}^{(j)}(f,T)$ is given explicitly by

$$\begin{eqnarray}{\mu }_{S}^{(j)}(f,T)\,=\, & \displaystyle \frac{1}{T}{\displaystyle \int }_{0}^{T}{\displaystyle \int }_{0}^{T}{{\rm{d}}{t}}_{1}{{\rm{d}}{t}}_{2}\cos (f({t}_{1}-{t}_{2})){\mathbb{E}}\{{X}_{{t}_{1}}^{(j)}{X}_{{t}_{2}}^{(j)}\},\\ & & \,\end{eqnarray} \tag{ 5 }$$

${\mathbb{E}}\{{X}_{{t}_{1}}^{(j)}{X}_{{t}_{2}}^{(j)}\}=2D\,\min ({t}_{1},{t}_{2})$ being the covariance of ${X}_{t}^{(j)}$. For BM, one then finds the standard result [1]

$$\begin{eqnarray}&&{\mu }_{S}^{(j)}(f,T)=\displaystyle \frac{4D}{{f}^{2}}\left[1-\displaystyle \frac{\sin ({fT})}{{fT}}\right],\end{eqnarray} \tag{ 6 }$$

such that

$$\begin{eqnarray}&&{\mu }_{S}^{(j)}(f,\infty )=\displaystyle \frac{4D}{{f}^{2}}.\end{eqnarray} \tag{ 7 }$$

Hence, for BM the standard PSD ${\mu }_{S}^{(j)}(f,\infty )$, defined as an average over an ensemble of trajectories, is described by a power-law with characteristic exponent β = 2 and an amplitude which is linearly proportional to the diffusion coefficient D.

Going beyond the text-book definition in equation (5), we now define the PSD of a single component ${X}_{t}^{(j)}$:

$$\begin{eqnarray}&&{S}_{T}^{(j)}(f)=\displaystyle \frac{1}{T}{\int }_{0}^{T}{\int }_{0}^{T}{{\rm{d}}{t}}_{1}{{\rm{d}}{t}}_{2}\cos (f({t}_{1}-{t}_{2})){X}_{{t}_{1}}^{(j)}{X}_{{t}_{2}}^{(j)},\end{eqnarray} \tag{ 8 }$$

and another property of interest here, the partial PSD of the trajectory ${{\bf{X}}}_{t}$

$$\begin{eqnarray}{\tilde{S}}_{T}^{(k)}(f) & = & \displaystyle \frac{1}{T}{\displaystyle \int }_{0}^{T}{\displaystyle \int }_{0}^{T}{{\rm{d}}{t}}_{1}{{\rm{d}}{t}}_{2}\cos (f({t}_{1}-{t}_{2}))\\ & & \times [{X}_{{t}_{1}}^{(1)}{X}_{{t}_{2}}^{(1)}+{X}_{{t}_{1}}^{(2)}{X}_{{t}_{2}}^{(2)}+\ldots +{X}_{{t}_{1}}^{(k)}{X}_{{t}_{2}}^{(k)}],\end{eqnarray} \tag{ 9 }$$

in which we take into account the contributions of k, (k = 1, 2, ..., d), components of a d-dimensional BM. Clearly, for k = 1 the definitions (8) and (9) coincide. The PSDs ${S}_{T}^{(j)}(f)$ and ${\tilde{S}}_{T}^{(k)}(f)$ are f- and T-parameterized random variables: the first moment ${\mu }_{S}(f,T)$ of a single-trajectory PSD (8) is given by the standard result (6), while the first moment of the partial PSD (9), due to statistical independence of the components, is given by expression (6) multiplied by k. Our goal is to evaluate exactly the full probability density function for ${\tilde{S}}_{T}^{(k)}(f)$.

To get an idea of how representative of the actual behavior of a single-trajectory PSD the result (6) is, we first look at the variance ${\sigma }_{S}^{2}(f,T)$ of a single-component single-trajectory PSD,

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{S}^{2}(f,T) & = & {\mathbb{E}}\{{({S}_{T}^{(j)}(f))}^{2}\}-{\mu }_{S}^{2}(f,T)\\ & = & \displaystyle \frac{20{D}^{2}}{{f}^{4}}\phantom{\rule{}{3.0ex}}[1-(6-\cos ({fT}))\displaystyle \frac{2\sin ({fT})}{5{fT}}+\displaystyle \frac{(17-\cos (2{fT})-16\cos ({fT}))}{10{f}^{2}{T}^{2}}\phantom{\rule{}{3.0ex}}]\\ & & \end{array}\end{eqnarray} \tag{ 10 }$$

Expression (10) permits us to determine the corresponding coefficient of variation

$$\begin{eqnarray}&&{\gamma }_{S}={\sigma }_{S}(f,T)/{\mu }_{S}(f,T)\end{eqnarray} \tag{ 11 }$$

of the yet unknown distribution of a single-component single-trajectory PSD.

In figure 1 we depict γ_S, which is a function of the product ${f}\,{T}$ exclusively. We observe that γ_S approaches $\sqrt{2}$ for fixed T and $f\to 0$, and tends to the asymptotic value $\sqrt{5}/2$ (thin horizontal line in figure 1) when $T\to \infty$ at any fixed f > 0. Overall, for any value of ${fT}$ the coefficient of variation appears to be greater than unity, meaning that the standard deviation ${\sigma }_{S}(f,T)$ is greater than the average, ${\mu }_{S}(f,T)$, which signals that the parental distribution is effectively 'broad', despite the fact that it evidently has well-defined moments of arbitrary order. As a consequence, the average described by equation (6) may indeed not be representative of the actual behavior of a single-trajectory PSD. This fully validates our quest for the distribution $P({\tilde{S}}_{T}^{(k)}(f))$.

Zoom In Zoom Out Reset image size

3.1. Brownian motion: moment-generating-function, distribution $P({\tilde{S}}_{T}^{(k)}(f))$ and its moments

Our first goal is to calculate the moment-generating function of ${\tilde{S}}_{T}^{(k)}(f)$ in equation (9), defined formally as the following Laplace transform

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }={{\rm{\Phi }}}_{\lambda }({\tilde{S}}_{T}^{(k)}(f))\equiv {\mathbb{E}}\{\exp (-\lambda \,{\tilde{S}}_{T}^{(k)}(f))\}.\end{eqnarray} \tag{ 12 }$$

To calculate Φ_λ exactly, it appears convenient to use Wiener's representation of a given Brownian path ${X}_{t}^{(j)}$ in the form of a Fourier series with random coefficients

$$\begin{eqnarray}&&{X}_{t}^{(j)}=\displaystyle \frac{\sqrt{2T}}{\pi }\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{{\zeta }_{n}^{(j)}}{(n-1/2)}\sin \left(\displaystyle \frac{(n-1/2)\pi t}{T}\right),\end{eqnarray} \tag{ 13 }$$

where ${\zeta }_{n}^{(j)}$ are independent, identically-distributed random variables with the distribution

$$\begin{eqnarray}&&P({\zeta }_{n}^{(j)})=\displaystyle \frac{1}{\sqrt{4\pi D}}\exp (-{({\zeta }_{n}^{(j)})}^{2}/4D).\end{eqnarray} \tag{ 14 }$$

The corresponding single-trajectory partial PSD ${\tilde{S}}_{T}^{(k)}(f)$ in equation (9) can be formally rewritten as

$$\begin{eqnarray}&&{\tilde{S}}_{T}^{(k)}(f)=\displaystyle \frac{2}{{\pi }^{2}}\displaystyle \sum _{j=1}^{k}\left[{\left(\displaystyle \sum _{n=1}^{\infty }{g}_{n}{\zeta }_{n}^{(j)}\right)}^{2}+{\left(\displaystyle \sum _{n=1}^{\infty }{h}_{n}{\zeta }_{n}^{(j)}\right)}^{2}\right],\end{eqnarray} \tag{ 15 }$$

where the functions ${g}_{n}={g}_{n}(f,T)$ and ${h}_{n}={h}_{n}(f,T)$ are given explicitly by

$$\begin{eqnarray}&&{g}_{n}=\displaystyle \frac{T}{(n-1/2)}\displaystyle \frac{{(-1)}^{n}{fT}\sin ({fT})+\pi (n-1/2)}{({\pi }^{2}{(n-1/2)}^{2}-{({fT})}^{2})},\end{eqnarray} \tag{ 16a }$$

and

$$\begin{eqnarray}&&{h}_{n}=-\displaystyle \frac{T}{(n-1/2)}\displaystyle \frac{{(-1)}^{n}{fT}\cos ({fT})}{({\pi }^{2}{(n-1/2)}^{2}-{({fT})}^{2})}.\end{eqnarray} \tag{ 16b }$$

Inserting expression (15) into (12) and using the identity

$$\begin{eqnarray}&&\exp \left(-\displaystyle \frac{{b}^{2}}{4a}\right)=\sqrt{\displaystyle \frac{a}{\pi }}{\int }_{-\infty }^{\infty }{\rm{d}}{x}\exp (-{{ax}}^{2}+{\rm{i}}{b}{x}),\end{eqnarray} \tag{ 17 }$$

we rewrite equation (12) in the factorized form

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }={\left[\displaystyle \frac{\pi }{8\lambda }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\rm{d}}{x}{\rm{d}}{y}{\rm{\exp }}\left(-\displaystyle \frac{{\pi }^{2}}{8\lambda }({x}^{2}+{y}^{2})\right){{\mathbb{E}}}_{\zeta }\left\{{\rm{\exp }}\left({\rm{i}}\displaystyle \sum _{n=1}^{\infty }({{xg}}_{n}+{{yh}}_{n}){\zeta }_{n}^{(j)}\right)\right\}\right]}^{k},\end{eqnarray} \tag{ 18 }$$

where the subscript ζ in the averaging operator ${{\mathbb{E}}}_{\zeta }\{\ldots \}$ signifies that averaging over all paths of the component ${X}_{t}^{(j)}$ is replaced by an equivalent operation, the averaging over all possible values of ${\zeta }_{n}^{(j)}$. Performing this averaging, as well as the integrations over dx and dy, we get

$$\begin{eqnarray}{{\rm{\Phi }}}_{\lambda }\,=\, & \left[1+\displaystyle \frac{8\lambda D}{{\pi }^{2}}\displaystyle \sum _{n=1}^{\infty }({g}_{n}^{2}+{h}_{n}^{2})+{\left(\displaystyle \frac{8\lambda D}{{\pi }^{2}}\right)}^{2}\right.{\left.\times \left(\displaystyle \sum _{n,l=1}^{\infty }{g}_{n}^{2}{h}_{l}^{2}-{\left(\displaystyle \sum _{n=1}^{\infty }{g}_{n}{h}_{n}\right)}^{2}\right)\right]}^{-k/2}.\\ & & \end{eqnarray} \tag{ 19 }$$

We focus next on the infinite sums entering equation (19). They can be calculated exactly, and expressed via the first and the second moments of a single-trajectory PSD,

$$\begin{eqnarray}&&\displaystyle \sum _{n=1}^{\infty }({g}_{n}^{2}+{h}_{n}^{2})\equiv \displaystyle \frac{{\pi }^{2}}{4D}{\mu }_{S}(f,T)\end{eqnarray} \tag{ 20a }$$

and

$$\begin{eqnarray}&&\displaystyle \sum _{n,l=1}^{\infty }{g}_{n}^{2}{h}_{l}^{2}-{\left(\displaystyle \sum _{n=1}^{\infty }{g}_{n}{h}_{n}\right)}^{2}\equiv {\left(\displaystyle \frac{{\pi }^{2}}{8D}\right)}^{2}{\mu }_{S}^{2}(f,T)(2-{\gamma }_{S}^{2}),\end{eqnarray} \tag{ 20b }$$

where the average and the variance (as well as the corresponding coefficient ${\gamma }_{S}$ of variation) of the single-component single-trajectory PSD are defined in equations (6), (10), and (11), respectively.

Consequently we realize that the moment-generating function in equation (19) can be cast into a more compact and physically meaningful form, which involves only the first moment and the variance (through the coefficient γ_S of variation) of the single-trajectory PSD,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }={[1+2{\mu }_{S}(f,T)\lambda +(2-{\gamma }_{S}^{2}){\mu }_{S}^{2}(f,T){\lambda }^{2}]}^{-k/2}.\end{eqnarray} \tag{ 21 }$$

This expression holds for any value of k, f and T.

Performing next the inverse Laplace transform of the function defined by equation (21), we arrive at the following expression for the desired probability density function of the single-trajectory partial PSD defined in equation (9), which also (as the result in equation (21)) holds for arbitrary f, arbitrary T and arbitrary number k of the projections of the trajectory ${{\bf{X}}}_{t}$ onto the coordinate axes,

$$\begin{eqnarray}P({\tilde{S}}_{T}^{(k)}(f)=S) & = & \displaystyle \frac{\sqrt{\pi }}{{2}^{\tfrac{k-1}{2}}{\rm{\Gamma }}(k/2)}\displaystyle \frac{{S}^{\tfrac{k-1}{2}}}{\sqrt{2-{\gamma }_{S}^{2}}{({\gamma }_{S}^{2}-1)}^{\tfrac{k-1}{4}}{\mu }_{S}^{\tfrac{k+1}{2}}(f,T)}\\ & & \times \exp \left(-\displaystyle \frac{1}{2-{\gamma }_{S}^{2}}\displaystyle \frac{S}{{\mu }_{S}(f,T)}\right){{I}}_{\tfrac{k-1}{2}}\left(\displaystyle \frac{\sqrt{{\gamma }_{S}^{2}-1}}{2-{\gamma }_{S}^{2}}\displaystyle \frac{S}{{\mu }_{S}(f,T)}\right).\end{eqnarray} \tag{ 22 }$$

Here, ${I}_{\alpha }(\ldots )$ is the modified Bessel function of the 1st kind and Γ(...) is the Gamma-function.

Before we proceed further, two remarks are in order. First, we note that the distribution (22) is the Bessel function distribution that has been used in mathematical statistics years ago as an example of a distribution with heavier than Gaussian tails (see, e.g., [75, 76]). Second, as already mentioned, for f = 0 the coefficient of variation is exactly equal to $\sqrt{2}$ and hence, the coefficient in front of the term quadratic in λ in expression (21) vanishes. In this particular case the distribution (22) simplifies to the χ²-distribution with k degrees of freedom, which is presented in appendix A.

Expression (22) permits us to straightforwardly calculate the moments of the partial single-trajectory PSD of arbitrary, not necessarily integer order $Q\gt -(k+1)/2$,

$$\begin{eqnarray}&&\displaystyle \frac{{\mathbb{E}}\{{({\tilde{S}}_{T}^{(k)}(f))}^{Q}\}}{{\mu }_{S}^{Q}(f,T)}=\displaystyle \frac{{\rm{\Gamma }}(Q+k){(2-{\gamma }_{S}^{2})}^{Q+k/2}}{{\rm{\Gamma }}(k)}\times {}_{2}{F}_{1}\left(\displaystyle \frac{Q+k}{2},\displaystyle \frac{Q+k+1}{2};\displaystyle \frac{k+1}{2};{\gamma }_{S}^{2}-1\right),\end{eqnarray} \tag{ 23 }$$

where ${}_{2}{F}_{1}(...)$ is the Gauss hypergeometric function. Note that the moments of order Q > 2 for an arbitrary k are all expressed through the first and the second (via γ_S) moments of the single-component single-trajectory PSD only, since the parental Gaussian process ${X}_{t}^{(j)}$ is entirely defined by its first two moments.

The distribution in equation (22) and the formula for the moments of order Q, equation (23), ensure that the single-trajectory PSD ${\tilde{S}}_{T}^{(k)}(f)$ has the form

$$\begin{eqnarray}&&{\tilde{S}}_{T}^{(k)}(f)={A}^{(k)}({\gamma }_{S})\,{\mu }_{S}(f,T),\end{eqnarray} \tag{ 24 }$$

where ${\mu }_{S}(f,T)$ is the first moment of this random variable, i.e., a deterministic function of f and T which sets the scale of variation of ${\tilde{S}}_{T}^{(k)}(f)$, and ${A}^{(k)}({\gamma }_{S})$ is a dimensionless random amplitude, whose moments of order Q are defined by the expression in the right-hand-side of equation (23). The relation in equation (24) has important conceptual consequences on which we will elaborate below.

3.2. Brownian motion: single-trajectory PSD

One infers from figure 1 that upon an increase of ${f}\,{T}$ the oscillatory terms in γ_S fade out, and γ_S saturates at the value ${\gamma }_{S}=\sqrt{5}/2$. Let us define ω_l as the value of ${f}\,{T}$ when the amplitude of the oscillatory terms in γ_S equals $\sqrt{5}/2+\varepsilon$, where ε > 0 is a small fixed number. Given that the amplitude of the oscillatory terms is a monotonically decreasing function of ${f}\,{T},$ one has that for ${f}\,{T}$ > ω_l the coefficient of variation $| {{\gamma }}_{{S}}({f}{T}\gt {{\omega }}_{{l}})-\sqrt{5}/2| \lt {\epsilon }$. Next, the decay law of the amplitude of oscillations can be readily derived from equations (6) and (10) to give that, in the leading in ε order, ${\omega }_{l}\sim 2/(5\varepsilon )$.

Then, for ${fT}\gt {\omega }_{l}$, up to terms proportional to ε, which can be made arbitrarily small, we see that the moments of the random amplitude ${A}^{(k)}({\gamma }_{S})$ in equation (24) are given by

$$\begin{eqnarray}&&{\mathbb{E}}\{{({A}^{(k)}({\gamma }_{S}))}^{Q}\}=\displaystyle \frac{{3}^{Q+k/2}{\rm{\Gamma }}(Q+k)}{{4}^{Q+k/2}{\rm{\Gamma }}(k)}\times {}_{2}{F}_{1}\left(\displaystyle \frac{Q+k}{2},\displaystyle \frac{Q+k+1}{2};\displaystyle \frac{k+1}{2};\displaystyle \frac{1}{4}\right),\end{eqnarray} \tag{ 25 }$$

meaning that in this limit ${A}^{(k)}({\gamma }_{S})$ becomes just a real number ${A}^{(k)}$ that is independent of f and T. This implies, in turn, that with any necessary accuracy, prescribed by the choice of ε, and for any realization of a trajectory ${{\bf{X}}}_{t}$,

$$\begin{eqnarray}&&{\tilde{S}}_{T}^{(k)}(f)=\displaystyle \frac{4D}{{f}^{2}}{A}^{(k)}(1+O(\varepsilon )),\,\mathrm{for}\ {f}\,{T}\in ({\omega }_{l},\infty ),\end{eqnarray} \tag{ 26 }$$

where the symbol O(ε) signifies that the omitted terms, (stemming from the oscillatory terms in ${\mu }_{S}(f,T)$, ${\sigma }_{S}^{2}(f,T)$ and hence, in ${\gamma }_{S}(f,T)$), have an amplitude smaller than ε.

Therefore, we arrive at the following conclusion, which is the main conceptual result of our analysis: for continuous-time BM and ${fT}\in ({\omega }_{l},\infty )$, the frequency-dependence of the PSD of any single trajectory is the same as that of the ensemble-averaged PSD with any desired accuracy, set by ε. Consequently, in response to the title question of our work, we conclude that for BM the 1/f²-dependence of the power spectrum can be deduced already from a single, sufficiently long trajectory, without any necessity to perform an additional averaging over an ensemble of such trajectories. In section 5 (see figure 7) below we show that this prediction made for BM—a somewhat idealized stochastic process with infinitesimal increments—holds indeed for discrete-time random walks and single-trajectory experiments, in which a BM trajectory is recorded at discrete instants of time and for a finite observation time.

Lastly, we note that since ${\mu }_{S}(f,\infty )$ is linearly proportional to the diffusion coefficient D and the latter appears to be multiplied by a random numerical amplitude, one cannot infer the ensemble-averaged diffusion coefficient from a single-trajectory PSD. The error in estimating D from a single trajectory is given precisely by the deviation of this amplitude from its averaged value. Below we discuss the fluctuations of this numerical amplitude.

3.3. Fluctuations of the amplitude ${A}^{(k)}$

The exact limiting probability density function of the amplitude ${A}^{(k)}$ in equation (26) follows from equation (22) and reads

$$\begin{eqnarray}&&P({A}^{(k)}=A)=\displaystyle \frac{2\sqrt{\pi }\,{A}^{\tfrac{k-1}{2}}}{\sqrt{3}\,{\rm{\Gamma }}(k/2)}\exp \left(-\displaystyle \frac{4}{3}A\right){{I}}_{\tfrac{k-1}{2}}\left(\displaystyle \frac{2}{3}A\right).\end{eqnarray} \tag{ 27 }$$

In figure 2 we depict this distribution for k = 1, 2 and 3 (curves from left to right). We observe that the very shape of the distribution depends on the number of components which are taken into account in order to evaluate the partial PSD. For k = 1 (that is, for BM in one-dimensional systems, or in two or three-dimensions but when only one of the components is being tracked), the distribution is a monotonically decreasing function with the maximal value at ${A}^{(1)}=0$. In contrast, for k = 2 and k = 3, $P({A}^{(k)})$ is a bell-shaped function with a left power-law tail and an exponential right tail.

Zoom In Zoom Out Reset image size

In figure 2 we also present a comparison of our analytical predictions in equation (27) against the results of numerical simulations. We use two methods to produce numerically the distributions of ${A}^{(k)}({\gamma }_{S})$ defined in equation (24) for the range of frequencies where ${f}^{2}{S}_{T}(f)$ is almost constant (see sections 4 and 5 for more details). The first method hinges on Wiener's representation of a BM in equation (13), which is truncated at the upper summation limit at N = 3.2 × 10⁶ (see section 4 below for more details). Crosses (×) in figure 2 depict the corresponding results obtained via averaging over 10⁵ trajectories generated using such a representation of BM. Further, we use a discrete-time representation of a BM—lattice random walks with unit spacing and stepping events at each tick of the clock—which are produced by Monte Carlos simulations. We set the observation time T = 1 and generate trajectories with $M={2}^{22}$ steps to get a single-trajectory periodogram. A thorough discussion of the domain of frequencies in which the periodogram yields the behavior specific to a continuous-time BM are presented below in section 5. Pluses (+) depict the results obtained numerically for random walks averaged over 10⁵ realizations of the process. Overall, we observe an excellent agreement between our analytical predictions, derived for BM—a continuous-time process with infinitesimal increments, and the numerical results. This signifies that the framework developed here is completely justified and can be used for a meaningful interpretation of stochastic trajectories obtained in single-trajectory experiments.

Moreover, the moments of the distribution in equation (27) are given by equation (25). We find that the average value is ${\mathbb{E}}\{{A}^{(k)}\}=k$, as it should be, and the variance $\mathrm{Var}({A}^{(k)})=5k/4$ so that the distribution broadens with increasing k. The coefficient of variation of the k-dependent distribution in equation (27) is equal to $\sqrt{5/(4k)}$ meaning that fluctuations become progressively less important for increasing k. We also remark that the most probable values for the cases k = 2 and k = 3 are ${A}_{{mp}}^{(2)}=3\mathrm{ln}(3)/4\approx 0.82$ and ${A}_{{mp}}^{(3)}\approx 1.74$, respectively, and are well below their average values.

Lastly, we quantify the effective broadness of the distribution in equation (27), focusing on its heterogeneity index

$$\begin{eqnarray}&&{\omega }^{(k)}=\displaystyle \frac{{A}_{1}^{(k)}}{{A}_{1}^{(k)}+{A}_{2}^{(k)}},\end{eqnarray} \tag{ 28 }$$

where ${A}_{1}^{(k)}$ and ${A}_{2}^{(k)}$ are the amplitudes drawn from two different independent realizations of ${{\bf{X}}}_{t}$. Such a diagnostic tool has been proposed in [77–79] in order to quantify fluctuations in the first passage phenomena in bounded domains, compare also with the discussion in [80, 81]. In our context, ${\omega }^{(k)}$ shows the likeliness of the event that two values of the amplitudes ${A}_{1}^{(k)}$ and ${A}_{2}^{(k)}$ will be equal to each other. With the distribution (27) of the random variable ${A}^{(k)}$ the distribution $P({\omega }^{(k)})$ of the heterogeneity index can be calculated via its integral representation [77, 78]

$$\begin{eqnarray}&&P({\omega }^{(k)}=\omega )=\displaystyle \frac{1}{{(1-\omega )}^{2}}{\displaystyle \int }_{0}^{\infty }{A}^{(k)}{{\rm{d}}{A}}^{(k)}\,P({A}^{(k)})\,P\left(\displaystyle \frac{\omega }{1-\omega }{A}^{(k)}\right).\end{eqnarray} \tag{ 29 }$$

Performing the integrals in relation (29), we find the following exact expressions for the distribution of the heterogeneity index: for k = 1 we have

$$\begin{eqnarray}P({\omega }^{(1)}=\omega ) & = & \displaystyle \frac{2\sqrt{3}}{\pi \sqrt{\omega (3+\omega (1-4\omega (2-\omega )))}}\\ & & \times \displaystyle \frac{{\bf{E}}(\exp (-\eta ))}{\sinh (\eta /2)},\,\eta =\mathrm{arcosh}\left(1+\displaystyle \frac{3}{2\omega (1-\omega )}\right),\end{eqnarray} \tag{ 30 }$$

where ${\bf{E}}(\ldots )$ is the complete elliptic integral of the second kind; for k = 2 the distribution has the form

$$\begin{eqnarray}&&P({\omega }^{(2)}=\omega )=\displaystyle \frac{2\omega (1-\omega )(39-20\omega (1-\omega ))}{{(3-2\omega )}^{2}{(1+2\omega )}^{2}},\end{eqnarray} \tag{ 31 }$$

while for k = 3 it follows that

$$\begin{eqnarray}&&P({\omega }^{(3)}=\omega )=-\displaystyle \frac{27}{16}\sqrt{\omega (1-\omega )}{\left[\displaystyle \frac{{d}^{3}}{{{dp}}^{3}}\exp \left(-\displaystyle \frac{3}{2}{\eta }^{* }\right){}_{2}{F}_{1}\left(\displaystyle \frac{1}{2},\displaystyle \frac{3}{2};2;\exp (-2{\eta }^{* })\right)\right]}_{p=1},\end{eqnarray} \tag{ 32 }$$

with

$$\begin{eqnarray}&&{\eta }^{* }=\mathrm{arcosh}\left(\displaystyle \frac{4{p}^{2}-{(1-\omega )}^{2}-{\omega }^{2}}{2\omega (1-\omega )}\right).\end{eqnarray} \tag{ 33 }$$

Note that the third-order derivative with respect to p in relation (32) can be taken explicitly producing, however, a rather cumbersome expression in terms of the complete elliptic integrals. For the sake of compactness, we nonetheless prefer the current notation.

The probability density functions in equations (30)–(32) are depicted in figure 3. We observe that for k = 1 the event in which two values of ${A}^{(1)}$ deduced from two different independent trajectories are equal to each other, (that is, when ${A}_{1}^{(1)}={A}_{2}^{(1)}$ and $\omega =1/2$), is the most unlikely, since it corresponds to the minimum of the distribution $P({\omega }^{(1)}=\omega )$ in equation (30). Therefore, for the case k = 1 the most probable outcome is that two values of ${A}^{(1)}$ obtained for two different realizations of BM will be very different from each other. For k = 2 and k = 3, $\omega =1/2$ corresponds to the most probable event but still the distributions in equations (31) and (32) appear to be rather broad so that a pronounced realization-to-realization variation of the amplitude is expected.

Zoom In Zoom Out Reset image size

In simulations of a continuous-time BM one often uses Wiener's representation (13), truncating it at the upper summation limit at some integer N,

$$\begin{eqnarray}&&{X}_{t}^{({\rm{tr}})}=\displaystyle \frac{\sqrt{2T}}{\pi }\displaystyle \sum _{n=1}^{N}\displaystyle \frac{{\zeta }_{n}}{(n-1/2)}\sin \left(\displaystyle \frac{(n-1/2)\pi t}{T}\right).\end{eqnarray} \tag{ 34 }$$

The sample paths of such a partial sum process are known to converge to the sample paths of the BM at the rate $1/\sqrt{N}$. Focusing on a single-component single-trajectory PSD in equation (8) (generalization to the d-dimensional case and a partial PSD ${\tilde{S}}_{T}^{(k)}(f)$ in equation (9) is straightforward) we have the following estimate for ${S}_{T}^{(j)}(f)$,

$$\begin{eqnarray}&&{S}_{T}^{({\rm{tr}})}(f)=\displaystyle \frac{2}{{\pi }^{2}}\left[{\left(\displaystyle \sum _{n=1}^{N}{g}_{n}{\zeta }_{n}\right)}^{2}+{\left(\displaystyle \sum _{n=1}^{N}{h}_{n}{\zeta }_{n}\right)}^{2}\right].\end{eqnarray} \tag{ 35 }$$

Below we examine how accurately ${S}_{T}^{({\rm{tr}})}(f)$ reproduces ${S}_{T}^{(j)}(f)$. We first consider the first moment of the truncated series and its variance, which are given by

$$\begin{eqnarray}&&{\mu }_{{\rm{tr}}}(f,T)=\displaystyle \frac{4D}{{\pi }^{2}}\displaystyle \sum _{n=1}^{N}({g}_{n}^{2}+{h}_{n}^{2})\end{eqnarray} \tag{ 36a }$$

and

$$\begin{eqnarray}&&{\sigma }_{{\rm{tr}}}^{2}(f,T)=\displaystyle \frac{32{D}^{2}}{{\pi }^{4}}\left[{\left(\displaystyle \sum _{n=1}^{N}{g}_{n}^{2}\right)}^{2}+{\left(\displaystyle \sum _{n=1}^{N}{h}_{n}^{2}\right)}^{2}\,+\,2{\left(\displaystyle \sum _{n=1}^{N}{g}_{n}{h}_{n}\right)}^{2}\right].\end{eqnarray} \tag{ 36b }$$

In figure 4 we present a comparison of the averaged single-component single-trajectory PSD and of its counterpart obtained from the truncated Wiener's series (34), as well as of the corresponding coefficients of variation of the two probability density functions. We observe a perfect agreement between the results obtained from the complete series and the truncated ones. Moreover, we see that ${\mu }_{{\rm{tr}}}(f,T)$ and ${\gamma }_{{\rm{tr}}}$ exhibit a uniform convergence to ${\mu }_{S}(f,T)$ and γ_S, respectively, for ${f}\,{T}\in (0,\pi N)$. This implies that keeping just N = 40 terms in the truncated Wiener's series permits us to describe reliably well the behavior of the pertinent properties over more than two decades of variation of f T. Extending this interval up to three decades will require keeping N ≈ 320 terms, for four decades N ≈ 3200 terms, and so on.

Zoom In Zoom Out Reset image size

We finally focus on the moment-generating function ${{\rm{\Phi }}}_{\lambda }({S}_{T}^{({\rm{tr}})}(f))$ of ${S}_{T}^{({\rm{tr}})}(f)$, and find

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }={[1+2{\mu }_{{\rm{tr}}}(f,T)\lambda +(2-{\gamma }_{{\rm{tr}}}^{2}){\mu }_{{\rm{tr}}}^{2}(f,T){\lambda }^{2}]}^{-1/2}.\end{eqnarray} \tag{ 37 }$$

We note that ${{\rm{\Phi }}}_{\lambda }({S}_{T}^{({\rm{tr}})}(f))$ has exactly the same form as the moment-generating function (21) (with k = 1) evaluated for a complete Wiener's series and hence, ${S}_{T}^{({\rm{tr}})}(f)$ has the distribution of exactly the same form as the one in equation (22) with the only difference that the first moment and the variance have to be replaced by their counterparts obtained via truncation of the Wiener's representation at some level N. Given that for ${fT}\in (0,\pi N)$ these properties are identical (see figure 4), it means that in this interval of variation of f T equation (37) coincides with (21), and the probability density function of ${S}_{T}^{({\rm{tr}})}(f)$ coincides with result (22). Outside of this interval, i.e., for ${f}\,{T}\gt \pi N$, the coefficient of variation of the truncated PSD jumps from $\sqrt{5}/2$ to $\sqrt{2}$ meaning that the distribution of ${S}_{T}^{({\rm{tr}})}(f)$ becomes the ${\chi }^{2}$-distribution (see appendix A), which is evidently a behavior.

In experiments or in Monte Carlo simulations of BM, the particle position is recorded at some discrete time moments such that one stores a given trajectory as a finite set of data. Here we analyze how the features unveiled in the previous sections will change, if instead of continuous-time BM we rather use a picture based on a discrete-time random walk.

Suppose that the time interval (0, T) is divided into M equally-sized subintervals Δ = T/M, such that the particle position is recorded at time moments ${t}_{m}={\rm{\Delta }}m$, $m=1,2,\,\ldots ,\,M$. We use the convention that at t = 0 the particle starts at the origin and focus on the behavior of a single-component PSD—the extension of our analysis over the general case of k components is straightforward but results in rather cumbersome expressions.

As a first step, we convert the integrals in equation (8) into the corresponding sums to get a periodogram

$$\begin{eqnarray}&&{R}_{M}(f)=\displaystyle \frac{{\bigtriangleup }^{2}}{T}\displaystyle \sum _{{m}_{\mathrm{1,2}}=0}^{M}\cos (f\bigtriangleup ({m}_{1}-{m}_{2}))\,{X}_{\bigtriangleup {m}_{1}}{X}_{\bigtriangleup {m}_{2}},\end{eqnarray} \tag{ 38 }$$

where we now denote a single-trajectory PSD as R_M(f) to emphasize that it is a different mathematical object as compared to the PSD in equation (8). Since we now have only a finite amount of points instead of a continuum, R_T(f) in equation (38) will be a periodic function of fT so that it may only approximate the behavior of the PSD for continuous-time BM in some range of frequencies at a given observation time. Further on, we use the scaling property of BM to rewrite the latter expression as

$$\begin{eqnarray}&&{R}_{M}(f)=\displaystyle \frac{2D{\bigtriangleup }^{3}}{T}\displaystyle \sum _{{m}_{\mathrm{1,2}}=0}^{M}\cos (f\bigtriangleup ({m}_{1}-{m}_{2}))\,{B}_{{m}_{1}}{B}_{{m}_{2}},\end{eqnarray} \tag{ 39 }$$

where B_m is now a trajectory of a standard lattice random walk with unit spacing and stepping at each clock tick m, m = 0, 1, 2, ..., M. Next, we write

$$\begin{eqnarray}&&{B}_{m}=\displaystyle \sum _{j=0}^{m}{s}_{j},\end{eqnarray} \tag{ 40 }$$

where s_j = ±1 are independent increments, and s₀ ≡ 0. Then, expression (39) becomes

$$\begin{eqnarray}&&{R}_{M}(f)=\displaystyle \frac{2D{{\rm{\Delta }}}^{2}}{M}\left[{\left(\displaystyle \sum _{j=1}^{M}{a}_{j}{s}_{j}\right)}^{2}+{\left(\displaystyle \sum _{j=1}^{M}{b}_{j}{s}_{j}\right)}^{2}\right],\end{eqnarray} \tag{ 41 }$$

where

$$\begin{eqnarray}{a}_{j} & = & \displaystyle \sum _{m=j}^{M}\cos (f{\rm{\Delta }}m)=\displaystyle \frac{1}{2}(\cos (f{\rm{\Delta }}j)+\cos (f{\rm{\Delta }}M)\\ & & +\cot (f{\rm{\Delta }}/2)(\sin (f{\rm{\Delta }}M)-\sin (f{\rm{\Delta }}j)))\end{eqnarray} \tag{ 42a }$$

and

$$\begin{eqnarray}{b}_{j} & = & \displaystyle \sum _{m=j}^{M}\sin (f{\rm{\Delta }}m)=\displaystyle \frac{1}{2}(\sin (f{\rm{\Delta }}j)+\sin (f{\rm{\Delta }}M)\\ & & +\cot (f{\rm{\Delta }}/2)(\cos (f{\rm{\Delta }}j)-\cos (f{\rm{\Delta }}M))).\end{eqnarray} \tag{ 42b }$$

Expression (41) is the discrete-time analog of the PSD in equation (15) obtained for the continuous-time BM.

5.1. Discrete-time case: mean and variance of a single-trajectory periodogram

At this point, it may be expedient to first look at the ensemble-averaged single-trajectory periodogram in equation (41) and at its variance, and to compare them against their continuous-time counterparts. Averaging the expression in equation (41) and its squared value over all possible realizations of the increments s_j, we have

$$\begin{eqnarray}&&{\mu }_{R}(f,T)=\displaystyle \frac{2D{{\rm{\Delta }}}^{2}}{M}\displaystyle \sum _{j=1}^{M}({a}_{j}^{2}+{b}_{j}^{2})=\displaystyle \frac{D{{\rm{\Delta }}}^{2}}{{\sin }^{2}(f{\rm{\Delta }}/2)}\left[1+\displaystyle \frac{1}{2\,M}-\displaystyle \frac{(\sin (f{\rm{\Delta }}M)+\sin (f{\rm{\Delta }}(M+1)))}{2\,M\,\sin (f{\rm{\Delta }})}\right],\end{eqnarray} \tag{ 43 }$$

and

$$\begin{eqnarray}{\sigma }_{R}^{2}(f,T) & = & \displaystyle \frac{8{D}^{2}{{\rm{\Delta }}}^{4}}{{M}^{2}}\left[{\left(\displaystyle \sum _{j=1}^{M}{a}_{j}^{2}\right)}^{2}-\displaystyle \sum _{j=1}^{M}{a}_{j}^{4}+{\left(\displaystyle \sum _{j=1}^{M}{b}_{j}^{2}\right)}^{2}-\displaystyle \sum _{j=1}^{M}{b}_{j}^{4}+2{\left(\displaystyle \sum _{j=1}^{M}{a}_{j}{b}_{j}\right)}^{2}-2\displaystyle \sum _{j=1}^{M}{a}_{j}^{2}\,{b}_{j}^{2}\right]\\ & = & \displaystyle \frac{20{D}^{2}{{\rm{\Delta }}}^{4}{\cos }^{4}(f{\rm{\Delta }}/2)}{{\sin }^{4}(f{\rm{\Delta }})}\left[1-\displaystyle \frac{1}{10M{\sin }^{2}(f{\rm{\Delta }})}(7-7\cos (2f{\rm{\Delta }})+6\cos (f{\rm{\Delta }}M)\right.\\ & & -\cos (2f{\rm{\Delta }}M)+\cos (2f{\rm{\Delta }}(M+1))\\ & & -6\cos (f{\rm{\Delta }}(M+2))+12\sin (f{\rm{\Delta }})\sin (f{\rm{\Delta }}M))\\ & & +\displaystyle \frac{1}{10{M}^{2}{\sin }^{2}(f{\rm{\Delta }})}(4+5\cos (2f{\rm{\Delta }})+8\cos (f{\rm{\Delta }})(1-\cos (f{\rm{\Delta }}M))\\ & & -2\cos (f{\rm{\Delta }}M)-2\cos (2f{\rm{\Delta }}M)+\cos (2f{\rm{\Delta }}(M+1))\\ & & -6\cos (f{\rm{\Delta }}(M+2))+16\sin (f{\rm{\Delta }})\sin (f{\rm{\Delta }}M))\phantom{\rule{}{3.5ex}}].\end{eqnarray} \tag{ 44 }$$

Note that these somewhat lengthy expressions (43) and (44) are exact, and valid for any M, Δ and f.

To illustrate the behavior of ${\mu }_{R}(f,T)$, ${\sigma }_{R}(f,T)$ and of the corresponding coefficient of variation ${\gamma }_{R}={\sigma }_{R}(f,T)/{\mu }_{R}(f,T)$ of the distribution of a single-trajectory periodogram R_M(f), we first depict in figure 5 these parameters as functions of fT for a rather small value of M, M = 10. We observe that they are periodic functions of fT with the prime period 2πM, (such that figure 1 presents a zoom of just one-half of the prime period) and exhibit peaks in the center, that is, at fT = πM, and at the end-points of the prime period, that is, for fT = 0 and fT = 2πM. More specifically, we see that in the middle of the prime period

$$\begin{eqnarray}&&{\mu }_{R}\left(\displaystyle \frac{\pi M}{T},T\right)=D{{\rm{\Delta }}}^{2}\end{eqnarray} \tag{ 45a }$$

and

$$\begin{eqnarray}&&{\sigma }_{R}^{2}\left(\displaystyle \frac{\pi M}{T},T\right)=2{D}^{2}{{\rm{\Delta }}}^{4}\left(1-\displaystyle \frac{2}{M}\right),\end{eqnarray} \tag{ 45b }$$

leading to

$$\begin{eqnarray}&&{\gamma }_{R}=\sqrt{2}\sqrt{1-\displaystyle \frac{2}{M}}.\end{eqnarray} \tag{ 45c }$$

For fixed T and $M\to \infty$, ${\mu }_{R}(\pi M/T,T)$ and ${\sigma }_{R}(\pi M/T,T)$ both tend to zero, but their ratio γ_R stays finite and tends to $\sqrt{2}$. Next, at the end-points of the prime period ${\mu }_{R}(f,T)$ and ${\sigma }_{R}(f,T)$ are given by

$$\begin{eqnarray}&&{\mu }_{R}\left(\displaystyle \frac{2\pi M}{T},T\right)=\displaystyle \frac{2{{DT}}^{2}}{3}\left(1+\displaystyle \frac{3}{2M}+\displaystyle \frac{1}{2{M}^{2}}\right)\end{eqnarray} \tag{ 46a }$$

and

$$\begin{eqnarray}&&{\sigma }_{R}^{2}\left(\displaystyle \frac{2\pi M}{T},T\right)=\displaystyle \frac{8{D}^{2}{T}^{4}}{9}\left(1+\displaystyle \frac{6}{5M}-\displaystyle \frac{5}{4{M}^{2}}+O\left(\displaystyle \frac{1}{{M}^{3}}\right)\right),\end{eqnarray} \tag{ 46b }$$

such that they both depend on the observation time T and diverge when $T\to \infty$. In this case γ_R obeys

$$\begin{eqnarray}&&{\gamma }_{R}=\sqrt{2}\left(1-\displaystyle \frac{9}{10M}+\displaystyle \frac{9}{200{M}^{2}}+O\left(\displaystyle \frac{1}{{M}^{3}}\right)\right),\end{eqnarray} \tag{ 46c }$$

and also stays finite as $M\to \infty$.

Zoom In Zoom Out Reset image size

Consider next the behavior of ${\mu }_{R}(f,T)$ and ${\sigma }_{R}(f,T)$ for fT = πM/2, that is, at one-quarter of the prime period, and also for fT = 2πM/3, that is, at one-third of the prime period. Here, we have

$$\begin{eqnarray}&&{\mu }_{R}\left(\displaystyle \frac{\pi M}{2T},T\right)=2D{{\rm{\Delta }}}^{2}\end{eqnarray} \tag{ 47a }$$

and

$$\begin{eqnarray}&&{\sigma }_{R}^{2}\left(\displaystyle \frac{\pi M}{2T},T\right)=5{D}^{2}{{\rm{\Delta }}}^{4}\left(1-\displaystyle \frac{12}{5M}\right),\end{eqnarray} \tag{ 47b }$$

producing

$$\begin{eqnarray}&&{\gamma }_{R}=\displaystyle \frac{\sqrt{5}}{2}\sqrt{1-\displaystyle \frac{12}{5M}},\end{eqnarray} \tag{ 47c }$$

and

$$\begin{eqnarray}&&{\mu }_{R}\left(\displaystyle \frac{2\pi M}{3T},T\right)=\displaystyle \frac{4D{{\rm{\Delta }}}^{2}}{3},\end{eqnarray} \tag{ 48a }$$

$$\begin{eqnarray}&&{\sigma }_{R}^{2}\left(\displaystyle \frac{2\pi M}{3T},T\right)=\displaystyle \frac{20{D}^{2}{{\rm{\Delta }}}^{4}}{9}\left(1-\displaystyle \frac{12}{5M}\right),\end{eqnarray} \tag{ 48b }$$

which imply that here, as well, γ_R is given by equation (47c). Therefore, we have that for the last two situations the first moment and the variance of a single-trajectory periodogram both tend to zero when ${\rm{\Delta }}\to 0$ (that is, when $M\to \infty$ at a fixed T), which is similar to the behavior at one-half of the prime period. In contrast, here γ_R tends to $\sqrt{5}/2$, which is the asymptotic coefficient of variation of the distribution of a single-trajectory PSD for the continuous-time BM.

In figure 6 we present a comparison between the averaged PSD (6) for the continuous-time BM and its discrete-time counterpart (43). We observe a nearly perfect agreement between these two representations for a wide range of variation of ${fT}$. We see that the periodicity of ${\mu }_{R}(f,T)$ starts to matter in a noticeable way only at ${fT}\approx 4\times {10}^{2}$ for M = 10³, at ${fT}\approx 3\times {10}^{3}$ for M = 10⁴ and only at fT ≈ 2 × 10⁴ for M = 10⁵, respectively. Below we will quantify the onset of the deviation between equations (43) and (6) more accurately.

Zoom In Zoom Out Reset image size

Next, we come back to figure 1 in which we compare the coefficient ${\gamma }_{S}$ of variation of the distribution of a single-trajectory PSD for the case of the continuous-time BM, and an analogous coefficient γ_R of variation of the distribution in the discrete-time case, calculated using our results (43) and (44) for M = 10⁴ and M = 10⁵. Here we observe that the agreement is even better and extends over a fairly large range of values of fT, until γ_R starts to show a strong oscillatory behavior close to the middle of the prime period (see also figure 5). Overall, it appears that γ_R can be set equal to $\sqrt{5}/2$ with an arbitrary accuracy dependent on the choice of ε (see section 3.2) on a bounded interval ${fT}\in ({\omega }_{l},\pi M-{\omega }_{l})$ and, by symmetry, for ${fT}\in (\pi M+{\omega }_{l},2\pi M-{\omega }_{l})$.

Therefore, in response to the title question of this work we conclude that one can indeed observe the behavior specific to the continuous-time BM (that is, the 1/f² dependence of the averaged PSD and, by virtue of equation (26), of a single-trajectory PSD) in the discrete-time settings. However, the spectrum must be analyzed on a finite interval (ω_l, ω_r) of variation of fT, which is bounded from below by ω_l (see section 3.2)—ensuring that γ_R becomes equal (with any necessary accuracy) to $\sqrt{5}/2$ and hence, does not contribute to the frequency-dependence of a single trajectory PSD—and by ω_r from above, when μ_R(f, T) starts to deviate from μ_S(f, T). Extending the arguments presented in section 3.2, we define ω_r as the value of fT at which ${\mu }_{R}(f,T)/{\mu }_{S}(f,T)=1+\varepsilon$. To calculate ω_r, we note that the ratio of two ensemble-averaged PSDs can be very accurately approximated by

$$\begin{eqnarray}&&\displaystyle \frac{{\mu }_{R}(f,T)}{{\mu }_{S}(f,T)}\approx \displaystyle \frac{{({fT})}^{2}}{4{M}^{2}{\sin }^{2}({fT}/2M)},\end{eqnarray} \tag{ 49 }$$

so that we find that ω_r = 2αM, where α obeys ${\alpha }^{2}/{\sin }^{2}(\alpha )=1+\varepsilon$. For small ε, the parameter α is given, to leading order in ε, by $\alpha =\sqrt{3\varepsilon }$. Recalling the definition of ω_l (see section 3.2), we thus see that the 1/f²-dependence can be observed once ${\omega }_{r}/{\omega }_{l}\gg 1$, which means that the number M of recorded positions of each trajectory has to obey $M\gg 1/(5\sqrt{3}{\varepsilon }^{3/2})$.

In figure 7 we plot the logarithm of a single-trajectory periodogram versus the logarithm of f T for five individual, randomly chosen trajectories of discrete-time lattice random walks and also five trajectories of a BM generated by the truncated Wiener's series (see section 4). For random walks, (which make a step of unit length at each tick of the clock, so that D = 1/2), we set the observation time T = 1 and have $M={2}^{22}\approx 4.2\times {10}^{6}$ discrete points within the unit time interval. For such a choice of parameters, we set the accuracy parameter ε = 0.01 to get ω_l = 40 and ${\omega }_{r}=1.45\times {10}^{6}$. We observe that, indeed, for all these five trajectories and $\mathrm{ln}f\in (\mathrm{ln}{\omega }_{l}=3.69,\mathrm{ln}{\omega }_{r}=14.19)$, (that is, for f varying over more than five decades), the relation ${\tilde{S}}_{T}^{(k)}(f)\sim 1/{f}^{2}$ holds. Next, for a BM generated by the truncated Wiener's representation we set T = 1, D = 0.01 and take N = 3.2 × 10⁶. Here, as well, for all five trajectories we observe a perfect agreement with our prediction in equation (26) for $\mathrm{ln}f\in (\mathrm{ln}{\omega }_{l}=3.69$, $\mathrm{ln}\pi N=14.98)$, that is, for f spanning over more than five decades.

Zoom In Zoom Out Reset image size

Lastly, we plot in figure 8 (upper panel) a logarithm of a single-trajectory periodogram versus the logarithm of f, for five trajectories of polystyrene beads in experiments performed within a flow cell (see below for more details). Each trajectory represents a one-dimensional projection of the bead motion, which is a BM with the diffusion coefficient $D\approx 0.365\,{\mu {\rm{m}}}^{2}\,{{\rm{s}}}^{-1}$, as estimated from the bead's squared displacement averaged over 150 trajectories. Note that the value of D inferred from the averaged PSD appears to be pretty close, $D\approx 0.373\,{\mu {\rm{m}}}^{2}\,{{\rm{s}}}^{-1}$. Again, we observe that the periodograms of the bead motion agree very well with our prediction in equation (26), and clearly follow the 1/f²-dependence for each individual trajectory. We stress again the similarly to the behavior of the MSD (2) [57, 58] to the amplitude scatter of a single-trajectory PSD as function of f: apart from some small fluctuations, it is remarkably constant.

Zoom In Zoom Out Reset image size

In the lower panel in figure 8, using a set of 150 individual experimental trajectories, we construct the distribution of the amplitude A in the 1/f²-dependence and compare it against our prediction in equation (27). The distributions of both one- and two-dimensional motions are shown, i.e., k = 1 (projection along a line) and k = 2 (projection along the imaging plane), respectively. Specifically, the distributions of the amplitudes are constructed by computing A = f²S(f)/(4 D) for 350 frequencies in the range 50 ≤ fT ≤ 400. We note that the agreement with equation (27) is quite impressive for both k = 1 and k = 2, with the size of the error bars being less than the size of the symbols.

In our experimental setup, we used a flow cell made up of a cover slip and a plastic slide. Two holes were drilled in the plastic slide to form the inlet and outlet of the chamber. Then the plastic slide and cover slip were attached using double-sided tape and sealed with nail polish, which was left to cure for 24 h. 1.2 μm polystyrene beads (SVP-10-5, Spherotech, Lake Forest, IL) were diluted in phosphate-buffered saline (PBS) with 0.05% Tween 20. The sample was agitated, centrifuged and resuspended in PBS with 1% bovine serum albumin (BSA) and 0.05% Tween 20. Subsequently the bead suspension was introduced in the flow cell chamber and the holes were sealed with nail polish, followed by immediate imaging. The beads were imaged in an inverted microscope equipped with a 40× objective (Olympus PlanApo, N.A. 0.95) and a sCMOS camera (Andor Zyla 4.2) operated at 100 frames per second. Each recorded sequence consisted of 4096 frames. Bead tracking in the plane was performed in LabView using a cross-correlation based tracking algorithm [82]. An image of a polystyrene bead is shown in the inset of figure 8. Several rings are observed in the image because the bead is intentionally not in focus, in order to increase tracking accuracy.

We close this subsection with the following remark. Estimating the diffusion coefficient or the temporal evolution of the MSD, one often slices a long trajectory into many smaller ones, which permits to create some statistical average. In this case, it seems to be appropriate for a stationary process because one here estimates either a constant (diffusion coefficient) or deals with a monotonically growing function of time. For the PSD the situation is much more subtle, because in the discrete-time case it is a periodic function of the parameter ω, with an oscillatory behavior within each prime period. Therefore, in this case one needs to fit the PSD only on a particular interval of variation of ω described below our equation (49). While slicing the trajectories permits one to create some statistical ensemble here, as well, this occurs at the expense of shrinking the interval [ω_l, ω_r], and it is not clear a priori if such a procedure can be beneficial or detrimental to the analysis. This question will be studied in detail elsewhere.

5.2. Discrete-time case: moment-generating function and distribution of a single-trajectory periodogram

The expression in equation (41) permits us to write formally the moment-generating function and the distribution of a single-trajectory periodogram as

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }({R}_{M}(f)=R)\equiv {\mathbb{E}}\{\exp (-\lambda {R}_{M}(f))\}=\displaystyle \frac{1}{{2}^{M}}\displaystyle \sum _{\{{s}_{j}\}}\exp \left(-\displaystyle \frac{2D{{\rm{\Delta }}}^{2}\lambda }{M}\left[{\left(\displaystyle \sum _{j=1}^{M}{a}_{j}{s}_{j}\right)}^{2}+{\left(\displaystyle \sum _{j=1}^{M}{b}_{j}{s}_{j}\right)}^{2}\right]\right)\end{eqnarray} \tag{ 50a }$$

and

$$\begin{eqnarray}&&P({R}_{M}(f)=R)=\displaystyle \frac{1}{{2}^{M}}\displaystyle \sum _{\{{s}_{j}\}}\delta \left(R-\displaystyle \frac{2D{{\rm{\Delta }}}^{2}}{M}\left[{\left(\displaystyle \sum _{j=1}^{M}{a}_{j}{s}_{j}\right)}^{2}+{\left(\displaystyle \sum _{j=1}^{M}{b}_{j}{s}_{j}\right)}^{2}\right]\right),\end{eqnarray} \tag{ 50b }$$

where the summation extends over all possible realizations of the sequence of the increments s_j. These expressions are exact, but they are of a little use, except for the (uninteresting) case when M is sufficiently small such that they can be written down in an explicit form by enumerations of all possible sequences of s_j.

We notice next that the weighted sums

$$\begin{eqnarray}&&{W}_{M}=\displaystyle \sum _{j=1}^{M}{a}_{j}{s}_{j}=\displaystyle \sum _{j=1}^{M}\left(\displaystyle \sum _{m=j}^{M}\cos (f{\rm{\Delta }}m)\right){s}_{j},\end{eqnarray} \tag{ 51a }$$

and

$$\begin{eqnarray}&&{V}_{M}=\displaystyle \sum _{j=1}^{M}{b}_{j}{s}_{j}=\displaystyle \sum _{j=1}^{M}\left(\displaystyle \sum _{m=j}^{M}\sin (f{\rm{\Delta }}m)\right){s}_{j},\end{eqnarray} \tag{ 51b }$$

appearing in equations (50a) and (50b) can be considered as two coupled discrete-time random walks with variable, generally irrational and incommensurate step-length. From this observation, some general conclusions can be reached about the expressions in equations (50a) and (50b): (i) for arbitrarily large M, both W_M and V_M are bounded from above by their maximal displacements, meaning that $P({R}_{M}(f)=R)$ has a bounded support and there is an upper cut-off value R_max(M) above which $P({R}_{M}(f))$ is identically zero. At the same time (apart from some special cases, for instance, when fΔ/π is a rational number, see appendices A and B) the trajectories W_M and V_M, for whatever large M, will never visit the origin again which signifies that $P({R}_{M}(f)=0)\equiv 0$. (ii) Moreover, we may expect that for an arbitrarily large M there is some gap R_min(M) below which $P({R}_{M}(f))$ is identically zero too, that is, W_M and V_M do not visit simultaneously some vicinity of the origin. These two specific features of the discrete-time settings make a significant difference as compared to the truncated Wiener's representation, in which analogous sums ${\sum }_{n=1}^{N}{g}_{n}{\zeta }_{n}$ and ${\sum }_{n=1}^{N}{h}_{n}{\zeta }_{n}$ are not bounded, and may, in principle, be equal to zero, since the increments ζ_n are continuous random variables with the support ($-\infty ,\infty$).

We proceed with the derivation of closed-form expressions for the moment-generating function and for the distribution of a single-trajectory periodogram. Using the integral identity in equation (17), we cast equations (50a) and (50b) into a form which permits to perform the averaging very directly. This leads to

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }({R}_{M}(f))=\displaystyle \frac{M}{8\pi D{{\rm{\Delta }}}^{2}\lambda }{\int }_{-\infty }^{\infty }{\rm{d}}{x}{\int }_{-\infty }^{\infty }{\rm{d}}{y}\exp \left(-\displaystyle \frac{M}{8D{{\rm{\Delta }}}^{2}\lambda }({x}^{2}+{y}^{2})\right)\displaystyle \prod _{j=1}^{M}\cos ({{xa}}_{j}+{{yb}}_{j})\end{eqnarray} \tag{ 52a }$$

and

$$\begin{eqnarray}&&P({R}_{M}(f)=R)=\displaystyle \frac{M}{8\pi D{{\rm{\Delta }}}^{2}}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{\rm{d}}{x}\,{\rm{d}}{y}\,{J}_{0}\left(\sqrt{\displaystyle \frac{M({x}^{2}+{y}^{2})R}{2D{{\rm{\Delta }}}^{2}}}\right)\displaystyle \prod _{j=1}^{M}\cos ({{xa}}_{j}+{{yb}}_{j}),\end{eqnarray} \tag{ 52b }$$

where J₀(...) is the Bessel function of the 1st kind. Expressions (52a) and (52b) are formally exact and hold for arbitrary M, f and Δ. We note that there is a set of magic frequencies f such that fΔ = πq with q a rational number, when a_j and b_j become periodic functions of j and the kernel ${\prod }_{j=1}^{M}\cos ({{xa}}_{j}+{{yb}}_{j})$ can be written down explicitly in form of a finite series of cosines. Then, the integrations in equations (52a) and (52b) can be performed exactly resulting in exact closed-form expressions for ${{\rm{\Phi }}}_{\lambda }({R}_{M}(f))$ and $P({R}_{M}(f)=R)$. Several examples of such calculations are presented in appendices A and B. However, for arbitrary f, analytical calculation of the integrals in results (52a) and (52b) is certainly beyond reach and one has to resort to some approximations.

5.3. Discrete-time case in the limit M ≫ 1: limiting forms of the moment-generating function and of the distribution of a single-trajectory periodogram

Meaningful approximations of expressions (52a) and (52b) can be obtained in the large-M limit. One has, however, to distinguish between the case when fT is fixed and $M\to \infty$, and the case when $M\to \infty$ but fT scales with M and may take any value within the interval (0, 2 πM).

5.3.1. Fixed fT and $M \rightarrow \infty$

For fixed T and not too large λ, the parameter $M/(8D{{\rm{\Delta }}}^{2}\lambda )={M}^{3}/(8{{DT}}^{2}\lambda )$ is large when M ≫ 1, so that the exponential function in equation (52a) is small everywhere except for the vicinity of the origin. Moreover, the integrals in equations (52a) and (52b) look very similar to the ones appearing in the theory of Pearson's random walks with a variable step length (see, e.g., [83]) and involve a kernel ${\prod }_{j=1}^{M}\cos ({{xa}}_{j}+{{yb}}_{j})$. It is well known that in the limit $M\to \infty$ the behavior of such a kernel is determined overwhelmingly by its behavior in the vicinity of x = 0 and y = 0, so that one has

$$\begin{eqnarray}&&\displaystyle \prod _{j=1}^{M}\cos ({{xa}}_{j}+{{yb}}_{j})\sim \exp \left(-\displaystyle \frac{1}{2}\displaystyle \sum _{j=1}^{M}{({{xa}}_{j}+{{yb}}_{j})}^{2}\right).\end{eqnarray} \tag{ 53 }$$

The accuracy of this approximation is checked in figure 9 in which we depict a density plot of the difference κ(x, y) of ${\prod }_{j=1}^{M}\cos ({{xa}}_{j}+{{yb}}_{j})$ and of the exponential function in the right-hand-side of (53). Already for M = 10² the difference κ(x, y) is zero almost everywhere and in the regions where it deviates from zero it is nonetheless numerically very small.

Zoom In Zoom Out Reset image size

Next, we have that for fT fixed and M ≫1 the ensemble-averaged periodogram and its variance admit the following representations

$$\begin{eqnarray}&&{\mu }_{R}(f,T)={\mu }_{S}(f,T)+\displaystyle \frac{4D{\sin }^{2}({fT}/2)}{{f}^{2}}\displaystyle \frac{1}{M}+O\left(\displaystyle \frac{1}{{M}^{2}}\right)\end{eqnarray} \tag{ 54a }$$

and

$$\begin{eqnarray}{\sigma }_{R}^{2}(f,T) & = & {\sigma }_{S}^{2}(f,T)+\displaystyle \frac{{D}^{2}}{{f}^{4}}\phantom{\rule{}{3.4ex}}(4(\cos (2{fT})-6\cos ({fT})-7)\\ & & \left.+8\left(7-\cos ({fT})\displaystyle \frac{\sin ({fT})}{{fT}}\right)\right)\displaystyle \frac{1}{M}+O\left(\displaystyle \frac{1}{{M}^{2}}\right),\end{eqnarray} \tag{ 54b }$$

where the omitted terms have a bounded amplitude for any value of fT and decay with the growth of M in proportion to the second inverse power of M.

Inserting expression (53) into equations (52a) and (52b) and performing the integrations, we find that in the limit M ≫ 1, for fixed f T and bounded λ, the moment-generating function and the distribution of a single trajectory periodogram are given up to terms of order $O(1/M)$, by our previous results (21) and (22).

5.3.2. Arbitrary fT ∈ (0, 2 πM) and $M \rightarrow \infty$

We now relax the condition that fT is fixed, and suppose that it can attain any value within the prime period 2πM. In other words, we consider the situation in which f T scales with M such that the expansions (54a) and (54b) are invalid, and we can no longer discard the correction terms. In this case, we find that the moment-generating function of a single-trajectory periodogram is given by

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }({R}_{M}(f))={\left[1+\displaystyle \frac{4D{{\rm{\Delta }}}^{2}\lambda }{M}\left(\displaystyle \sum _{j=1}^{M}({a}_{j}^{2}+{b}_{j}^{2})\right)\lambda +\displaystyle \frac{16{D}^{2}{{\rm{\Delta }}}^{4}}{{M}^{2}}\left(\left(\displaystyle \sum _{j,i=1}^{M}{a}_{j}^{2}{b}_{i}^{2}\right)-{\left(\displaystyle \sum _{j=1}^{M}{a}_{j}{b}_{j}\right)}^{2}\right){\lambda }^{2}\right]}^{-1/2}.\end{eqnarray} \tag{ 55 }$$

Differentiating this expression once and twice, and then setting λ = 0, we observe that (55) correctly reproduces the first moment of the periodogram, but yields an incorrect expression for the variance, meaning that in this limit the approximation (53) is insufficient (although it works fairly well in the previous case with fT kept fixed). In other words, in order to obtain a moment-generating function which reproduces correctly first two moments of a single-trajectory periodogram, one has to go beyond this approximation. For Pearson's random walks with j-dependent step-lengths a_j and b_j this turns out to be quite a complicated problem which we are not in the position to solve here.

Conversely, on intuitive grounds, we may conjecture such a form of the moment-generating function which reproduces the first and the second moment correctly. For $M\to \infty$, this is given by

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }={[1+2{\mu }_{R}(f,T)\lambda +(2-{\gamma }_{R}^{2}){\mu }_{R}^{2}(f,T){\lambda }^{2}]}^{-1/2},\end{eqnarray} \tag{ 56 }$$

that is, has precisely the same form as the moment-generating function in the continuous-time case, equation (21), but with the first two moments replaced by the analogous properties of a single-trajectory periodogram. Evidently, a generalization to the case of a k-component periodogram amounts to merely replacing 1/2 by k/2.

In turn, the form (56) implies that the distribution of a single-trajectory periodogram is given by equation (22) with ${\mu }_{S}(f,T)$ replaced by ${\mu }_{R}(f,T)$ and ${\sigma }_{S}^{2}(f,T)$ replaced by its discrete-time counterpart ${\sigma }_{R}^{2}(f,T)$, that is,

$$\begin{eqnarray}&&P({R}_{M}(f)=R)=\displaystyle \frac{1}{\sqrt{2-{\gamma }_{R}^{2}}{\mu }_{R}(f,T)}\exp \left(-\displaystyle \frac{1}{2-{\gamma }_{R}^{2}}\displaystyle \frac{R}{{\mu }_{R}(f,T)}\right){{I}}_{0}\left(\displaystyle \frac{\sqrt{{\gamma }_{R}^{2}-1}}{2-{\gamma }_{R}^{2}}\displaystyle \frac{R}{{\mu }_{R}(f,T)}\right).\end{eqnarray} \tag{ 57 }$$

We are unable to prove expression (56) (and hence, equation (57)) in the general case of arbitrary f. However, as we have already remarked, there are plenty of cases in which the moment-generation function and the distribution in equations (52a) and (52b) can be calculated exactly. In appendix B we present such calculations for several particular values of the frequency, such that $f{\rm{\Delta }}=2\pi /3$, fΔ = π/2, fΔ = 2 π and fΔ = π, and show that the exact discrete-time results do indeed converge to the asymptotic forms (56) and (57). Therefore, we have all grounds to believe that expression (57) is correct.

Lastly, we note that equation (57) permits us to make a substantial generalization of the single-trajectory relation (26). We see that for ${fT}\in ({\omega }_{l},\pi M-{\omega }_{l})$ (in which domain ${\gamma }_{R}=\sqrt{5}/2$ with any necessary accuracy set by the choice of ε), a single-trajectory periodogram obeys

$$\begin{eqnarray}&&{R}_{M}(f)=A\,{\mu }_{R}(f,T),\end{eqnarray} \tag{ 58 }$$

implying that the spectrum of a single-trajectory periodogram should be the same as for the ensemble-averaged periodogram, while only the amplitude A will be a fluctuating property. The limiting distribution of the amplitude is given by equation (27).

We studied here in detail the statistical properties of the PSD of a single trajectory of a d-dimensional BM. We calculated exactly the moment-generating function, the full probability density function and the moments of arbitrary, not necessarily integer order of such a PSD in the most general case of arbitrary frequency f, arbitrary (not necessarily infinite) observation time T and arbitrary number of projections of a given trajectory onto the coordinate axes. We showed that for a sufficiently large T (and the frequency f > 0) a single-trajectory PSD for any realization of the process is proportional to its first moment, which embodies the full dependence on the frequency and on the diffusion coefficient. This implies that the correct frequency-dependence specific to an ensemble of trajectories can be deduced already from a single trajectory and solely the numerical proportionality factor in this relation is random, and varies from realization to realization. Due to this fact, one cannot infer the precise value of the ensemble-averaged diffusion coefficient from a single-trajectory PSD since the ensemble-averaged PSD is linearly proportional to the diffusion coefficient which thus appears to be multiplied by a random amplitude. The distribution function of this amplitude was also calculated exactly here and its effective width was quantified using standard criteria.

Moreover, we addressed several questions emerging in connection with the numerical and experimental verification of our analytical predictions for the continuous-time BM. To this end, we first considered Wiener's representation of BM in form of an infinite Fourier series with random coefficients, whose truncated version is often used in numerical simulations. We showed that the distribution of the single-trajectory PSD obtained from Wiener's series in which just N terms are kept, instead of an infinite number, has exactly the same form as the one obtained for the continuous-time BM when fT is within the interval (0, πN). Outside this interval, the probability density function of the truncated PSD converges to a different form.

Next, we examined the case when a trajectory of continuous-time BM is recorded at some discrete time moments, so that the whole trajectory is represented by a set of M points and the PSD, called a periodogram, becomes a periodic function of the product fT with the prime period equal to 2πM. We analyzed several aspects of this discrete-time problem. Namely, we studied how big M should be taken at a fixed observation time T so that we may recover the results obtained for the continuous-time BM. Apart of that, we studied the limiting forms of the distribution of a single-trajectory periodogram and showed, in particular, that for fT kept fixed and $M\to \infty$, the latter converges to the form obtained for the continuous-time BM. In contrast, when fT is left arbitrary so that it may assume any value within the prime period, that is, fT ∈ (0, 2πM), the limiting distribution of a single-trajectory periodogram converges to a different form as $M\to \infty$. Our analysis revealed the remarkable observation that for fT belonging to a certain interval within the prime period, a single-trajectory periodogram equals, up to a random numerical amplitude, the ensemble-averaged periodogram, and the latter embodies the full dependence on f and T. Therefore, similarly to the continuous-time case, the correct spectrum can be obtained already from a single trajectory.

To check our analytical predictions we performed numerical analyzes, using a truncated Wiener's representation and also Monte Carlo simulations of discrete-time random walks, as well as an experimental analysis of BM of polystyrene beads in aqueous solution. We confirmed our prediction of the relation connecting a single-trajectory PSD and its ensemble averaged counter-part and demonstrated that the former indeed exhibits a well-defined 1/f²-dependence, specific to a standard, ensemble-averaged PSD. Furthermore, our numerical analysis confirmed the form of the distribution function of the amplitude in this relation.

The theoretical analysis for the case of BM presented here can be extended in several directions. In particular, one may inquire about analogous distributions for anomalous diffusion processes [57], as exemplified, for instance, by superdiffusive Lévy motion or subdiffusive continuous-time random walks with a broad distribution of the waiting times. Another important example of transport processes is given by a wide and experimentally relevant class of anomalous diffusion processes called fractional BM, for which the generalization of our analysis is relatively straightforward. Subdiffusive fractional BM is related to the overdamped, generalized Langevin equation with a power-law memory kernel typical for viscoelastic systems [84]. However, superdiffusive fractional BM was identified for active transport in biological cells [85]. We also mention the very active field of diffusion with varying diffusion coefficients. These can be systematically varying in space [64], time [86, 87], or be randomly varying in time [88–90] or space [91–93]. Lastly, a challenging field of further research is to search for such periodic functions of the frequency in equation (1), in place of the exponential function, for which a single-trajectory PSD will possess an ergodic property so that its variance will tend to zero in the limit $T\to \infty$.

DK, RM and GO wish to thank for warm hospitality the International Center for Advanced Studies in Buenos Aires, where this work was initiated. DK acknowledges a partial support from the NSF Grant No. 1401432. EM acknowledges funding from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement No. 694925).

We focus in this appendix on the moment-generating function and on the probability density function of a single-trajectory PSD for the continuous-time BM and of the discrete-time single-trajectory periodogram in the special case f = 0.

Consider first the case of BM. The first moment and the variance of the PSD for f = 0 follow from equations (6) and (10) and are given explicitly by

$$\begin{eqnarray}&&{\mu }_{S}(0,T)=\displaystyle \frac{2{{DT}}^{2}}{3},\end{eqnarray} \tag{ A1 }$$

and

$$\begin{eqnarray}&&{\sigma }_{S}^{2}(0,T)=\displaystyle \frac{8{D}^{2}{T}^{4}}{9},\end{eqnarray} \tag{ A2 }$$

so that the coefficient of variation is ${\gamma }_{S}\equiv \sqrt{2}$. Hence, in this case the coefficient in front of the quadratic term in equation (21) vanishes and the moment-generating function obeys

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }({\tilde{S}}_{T}^{(k)}(0))={\left[1+\displaystyle \frac{4{{DT}}^{2}}{3}\lambda \right]}^{-k/2}.\end{eqnarray} \tag{ A3 }$$

Taking the inverse Laplace transform, we then find that in this case the distribution $P({\tilde{S}}_{T}^{(k)}(0))$ is given by

$$\begin{eqnarray}&&P({\tilde{S}}_{T}^{(k)}(0)=S)={\left(\displaystyle \frac{3}{4{{DT}}^{2}}\right)}^{k/2}\displaystyle \frac{{S}^{k/2-1}}{{\rm{\Gamma }}(k/2)}\exp \left(-\displaystyle \frac{3S}{4{{DT}}^{2}}\right),\end{eqnarray} \tag{ A4 }$$

that is, it is the χ²-distribution with k degrees of freedom. Note that in contrast to the distribution (22) with f bounded away from zero, the one in equation (A4) does not attain a limiting form as the observation time $T\to \infty$. For any fixed S > 0 one has

$$\begin{eqnarray}&&P({\tilde{S}}_{T}^{(k)}(0)=S)\sim \displaystyle \frac{1}{{T}^{k}},\,\mathrm{as}\ T\to \infty .\end{eqnarray} \tag{ A5 }$$

In the discrete-time case, for f = 0, the first moment ${\mu }_{R}(0,T)$ of the periodogram, the variance ${\sigma }_{R}^{2}(0,T)$ and the variation coefficient γ_R are given by equations (46a)–(46c), respectively. We notice that they converge, as $M\to \infty$, to their continuous-time counterparts in equations (A1) and (A2), and ${\gamma }_{R}\to {\gamma }_{S}\equiv \sqrt{2}$, meaning that in this limit the moment-generating function and the distribution of the periodogram converge to the forms in equations (A3) and (A4).

For f = 0 in the discrete-time case, the moment-generating function and the distribution of the periodogram can be calculated exactly for arbitrary M, which permits us to estimate the rate of convergence. We constrain our analysis to the case when M is divisible by 4. We impose such a constraint just for a simplicity of exposition, exact solutions can be also obtained for other values of M resulting, however, in rather cumbersome expressions without providing any additional insight.

We notice that for f = 0 the coefficients b_j in equation (42b) are identically equal to zero, which implies that V_M ≡ 0, (51b), while the coefficients a_j in equation (42a) are explicitly given by ${a}_{j}\equiv M+1-j$. Consequently, the sum ${W}_{M}={\sum }_{j=1}(M+1-j){s}_{j}$ appearing in equations (50a) and (50b) can be thought of as a 'random walk' with a linearly 'shrinking' time-dependent step length (see, e.g., [94, 95] for other examples of such random walks). In this case, the product of cosines in equations (52a) and (52b) can be explicitly written down as

$$\begin{eqnarray}\displaystyle \prod _{j=1}^{M}\cos ({a}_{j}x+{b}_{j}y) & = & \displaystyle \prod _{j=1}^{M}\cos ((M+1-j)x)\\ & = & \displaystyle \frac{1}{{2}^{M}}\displaystyle \sum _{m=-M(M+1)/4}^{M(M+1)/4}{q}_{M}\left(\displaystyle \frac{M(M+1)}{4}-m\right)\cos (2{mx}),\end{eqnarray} \tag{ A6 }$$

where the coefficient q_M(j) is defined by

$$\begin{eqnarray}&&{q}_{M}(j)=\displaystyle \frac{1}{j!}{\left.\displaystyle \frac{{{\rm{d}}}^{j}}{{{\rm{d}}{q}}^{j}}{(-q,q)}_{M}\right|}_{q=0},\end{eqnarray} \tag{ A7 }$$

with ${(-q,q)}_{M}={\prod }_{m=1}^{M}(1+{q}^{m})$ being the q-Pochhammer symbol [96]. In other terms, q_M(j) is the numerical coefficient before the term q^j in the polynomial ${\prod }_{m=1}^{M}(1+{q}^{m})$. For $M=\infty$, ${q}_{\infty }(j)$ is simply the number of all possible partitions of an integer j into distinct parts. Note also that q_M(j) is symmetric around m = 0 such that ${q}_{M}(M(M+1)/4-m)={q}_{M}(M(M+1)/4+m)$, $m=1,2,\,\ldots ,\,M(M+1)/4$.

Inserting next expression (A6) into equations (52a) and (52b) and performing the integrations, we arrive at the following exact results for the moment-generating function and the distribution of a single-trajectory periodogram R_M(f) for f = 0:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }=\displaystyle \frac{1}{{2}^{M}}\displaystyle \sum _{m=-M(M+1)/4}^{M(M+1)/4}{q}_{M}\left(\displaystyle \frac{M(M+1)}{4}-m\right)\exp \left(-\displaystyle \frac{8D{{\rm{\Delta }}}^{2}{m}^{2}}{M}\lambda \right),\end{eqnarray} \tag{ A8 }$$

and

$$\begin{eqnarray}P({R}_{M}(0)=R) & = & \displaystyle \frac{1}{{2}^{M}}\left[{q}_{M}\left(\displaystyle \frac{M(M+1)}{4}\right)\delta (R)\right.\\ & & \left.+2\displaystyle \sum _{m=1}^{M(M+1)/4}{q}_{M}\left(\displaystyle \frac{M(M+1)}{4}-m\right)\delta \left(R-\displaystyle \frac{8D{{\rm{\Delta }}}^{2}{m}^{2}}{M}\right)\right].\end{eqnarray} \tag{ A9 }$$

Note that R has a discrete support on the interval $(0,{R}_{\max })$ with ${R}_{\max }=D{{\rm{\Delta }}}^{2}M{(M+1)}^{2}/2$. Note, as well, that the presence of a delta-function at R = 0 is the consequence of the choice of M; for M not divisible by 4 the coefficient in front of δ(R) equals zero because the 'random walk' W_M never visits the origin.

In figure A1, upper panel, we present a comparison of the expressions for the moment-generating functions in equation (A3) for the continuous-time case and in equation (A8) for the discrete-time case, for different densities of the discrete-time points for different values of T. We observe that for T = 1 and M = 12 (corresponding to Δ = 1/12) the agreement between the continuous- and the discrete-time results is fairly good up to λ = 10³. For T = 10 and M = 36, such that on average 3.6 points are recorded for a unit of time, the discrete-time result agrees with the continuous-time one up to $\lambda =3\times {10}^{2}$. Increasing the number of points up to M = 60, so that we have now 6 (instead of 3.6) points per unit of time, a perfect agreement between equations (A3) and (A8) extends up to λ = 10³. Lastly, for the observation time T = 20 with 60 recorded points (which corresponds to just 3 points per unit of time T) we have a perfect agreement between the discrete- and continuous-time predictions for values of λ up to λ = 3 × 10² and a significant departure for λ exceeding these values. In contrast, for T = 20 and M = 10², such that we have 5 points for each unit of time, a perfect agreement is observed for the whole range of variation of λ. We thus conclude that, in fact, once we seek to achieve a good agreement between the continuous- and discrete-time expressions for the moment-generation function, we do not need a very high precision in approximating the continuous-time trajectory by a discrete-time one.

Zoom In Zoom Out Reset image size

Next, in the lower panel of figure A1 we compare the probability density functions obtained for the continuous-time and the discrete-time cases. Since the distribution in the discrete-time case is a finite sum of delta-functions (see equation (A9)) it is appropriate to compare not the distributions themselves, but rather their integrated forms, the cumulative distribution functions ${ \mathcal F }$, obtained by integrating the forms in equations (A4) and (A9) over dR from 0 to R. From figure A1 (lower panel) we conclude that here too just 16 recorded points per unit of time appear to be enough to get a fairly good agreement between the distributions in the continuous- and the discrete-time cases.

In this appendix we present several stray examples for which the moment-generating function in equation (50a) and the distribution of a single-trajectory periodogram in equation (50b) can be calculated exactly. Our aim here is to demonstrate that in all these particular cases the exact moment-generating function and the distribution converge to the forms (56) and (57) as $M\to \infty$, which validates our conjecture. These particular choices of the frequency within the prime period are: fΔ = π, fΔ = π/2, fΔ = 2 π/3 and the endpoint of the prime period, fΔ = 2 π. We note that for all choices of the frequency the first moment μ_R(f, T) and the variance σ_R(f, T) converge as $M\to \infty$ to different values, as compared to their continuous-time counterparts. This means, in turn, that for fT within the prime period these are the forms in equations (56) and (57) which describe the moment-generating function and the distribution of the periodogram in the limit $M\to \infty$, but not the continuous-time results derived in section 3.

B.1. The case fΔ = π

We start with the simplest case, where the convergence of the moment-generating function (50a) and of the distribution of a single-trajectory periodogram (50b) to the limiting forms in equations (56) and (57) can be established analytically. We will again assume, for simplicity, that M is divisible by 4. We hasten to remark that such a constraint is not crucial—an exact solution can be obtained in the general case of an arbitrary M, but the resulting expressions will be more cumbersome.

For fΔ = π, the coefficients b_j in equation (42b) all vanish, such that the 'random walk' V_M (51b) stays at the origin. The coefficients a_j in equation (42a) are periodic with period 2; that being, ${a}_{j}={a}_{j+2}$, and are given by ${a}_{j}=(1+{(-1)}^{j})/2$ when M is divisible by 4. Consequently, the random walk W_M in equation (51a) steps at even time moments and pauses at the odd ones. In this case, the first moment and the variance of a single-trajectory periodogram are, respectively, given by equations (45a) and (45b), while the coefficient of variation obeys equation (45c) and tends to $\sqrt{2}$ as $M\to \infty$. This means that the limiting form of the distribution of a single trajectory periodogram for such a choice of f should be the χ²-distribution. We show below that this is indeed the case.

For M divisible by 4, the kernel in (52a) is given by

$$\begin{eqnarray}&&\displaystyle \prod _{j=1}^{M}\cos ({a}_{j}x+{b}_{j}y)={\cos }^{\tfrac{M}{2}}(x)=\displaystyle \frac{1}{{2}^{\tfrac{M}{2}}}\displaystyle \sum _{m=-\tfrac{M}{4}}^{\tfrac{M}{4}}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{2}}{\tfrac{M}{4}+m}\right)\,\cos (2\,{mx}),\end{eqnarray} \tag{ B1 }$$

where $\left(\genfrac{}{}{0em}{}{a}{b}\right)$ is the binomial coefficient. Substituting equation (B1) into equation (52a) and integrating term by term, we find that for fΔ = π the moment-generating function is given by the following sum of exponential functions:

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }\left({R}_{M}\left(\displaystyle \frac{\pi }{{\rm{\Delta }}}\right)\right)=\displaystyle \frac{1}{{2}^{\tfrac{M}{2}}}\displaystyle \sum _{m=-\tfrac{M}{4}}^{\tfrac{M}{4}}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{2}}{\tfrac{M}{4}+m}\right)\,\exp \left(-\displaystyle \frac{8D{{\rm{\Delta }}}^{2}{m}^{2}}{M}\lambda \right),\end{eqnarray} \tag{ B2 }$$

and hence, the distribution of a single-trajectory periodogram is given explicitly as a sum of delta-functions

$$\begin{eqnarray}&&P\left({R}_{M}\left(\displaystyle \frac{\pi }{{\rm{\Delta }}}\right)=R\right)=\displaystyle \frac{1}{{2}^{\tfrac{M}{2}}}\left[\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{2}}{\tfrac{M}{4}}\right)\delta (R)\,+\,2\displaystyle \sum _{m=1}^{\tfrac{M}{4}}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{2}}{\tfrac{M}{4}+m}\right)\delta \left(R-\displaystyle \frac{8D{{\rm{\Delta }}}^{2}}{M}{m}^{2}\right)\right].\end{eqnarray} \tag{ B3 }$$

Note that again the delta-peak at R = 0 is a consequence of the specific choice of M. For M not divisible by 4 such a peak is absent, precisely as it happens for a standard random walk, for which the parity of the number M of steps matters.

In the limit $M\to \infty$ it is legitimate to replace the summation operation by an integration and use

$$\begin{eqnarray}&&{2}^{-\tfrac{M}{2}}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{2}}{m}\right)\to \displaystyle \frac{2}{\sqrt{\pi M}}{{\rm{e}}}^{-\tfrac{{(m-M/4)}^{2}}{M/4}},\end{eqnarray} \tag{ B4 }$$

to get from expression (B2) the following limiting expression for the moment-generating function

$$\begin{eqnarray}{{\rm{\Phi }}}_{\lambda }\left({R}_{M}\left(\displaystyle \frac{\pi }{{\rm{\Delta }}}\right)\right) & \sim & {\displaystyle \int }_{-M/4}^{M/4}{\rm{d}}m\,\displaystyle \frac{2}{\sqrt{\pi M}}{{\rm{e}}}^{-(1+2D{{\rm{\Delta }}}^{2}\lambda )\tfrac{4{m}^{2}}{M}}\\ & = & \displaystyle \frac{{\rm{erf}}(\sqrt{M(1+2D{{\rm{\Delta }}}^{2}\lambda )}/2)}{\sqrt{1+2D{{\rm{\Delta }}}^{2}\lambda }}.\end{eqnarray} \tag{ B5 }$$

Recalling that ${\mu }_{R}(\pi /{\rm{\Delta }},T)=D{{\rm{\Delta }}}^{2}$ we see then that for $M\to \infty$ the right-hand-side of the latter equation converges to

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }\left({R}_{M}\left(\displaystyle \frac{\pi }{{\rm{\Delta }}}\right)\right)\sim {\left(1+2{\mu }_{R}\left(\displaystyle \frac{\pi }{{\rm{\Delta }}},T\right)\lambda \right)}^{-1/2},\end{eqnarray} \tag{ B6 }$$

which is precisely our conjectured expression (56) corresponding to ${\gamma }_{R}=\sqrt{2}$. The distribution is then obtained upon inverting the Laplace transform to get

$$\begin{eqnarray}&&P\left({R}_{M}\left(\displaystyle \frac{\pi }{{\rm{\Delta }}}\right)=R\right)=\displaystyle \frac{1}{\sqrt{2\pi {\mu }_{R}\left(\tfrac{\pi }{{\rm{\Delta }}},T\right)\,R}}\exp \left(-\displaystyle \frac{R}{2{\mu }_{R}\left(\tfrac{\pi }{{\rm{\Delta }}},T\right)}\right),\end{eqnarray} \tag{ B7 }$$

which is the χ²-distribution conjectured in equation (57). The convergence of equation (B2)–(B6), and also of (B3)–(B7) for different values of M is demonstrated in figure B1.

Zoom In Zoom Out Reset image size

B.2. The case $f{\rm{\Delta }}=\pi /2$

Consider the case when fT is at one-quarter of the prime period. In this case the coefficients a_j and b_j are periodic functions of j with period 4, i.e., ${a}_{j}={a}_{j+4}$ and ${b}_{j}={b}_{j+4}$. Within the period, their values are given by

$$\begin{eqnarray}\begin{array}{ccccccc} & & & & & & \\ & j & 1 & 2 & 3 & 4 & \\ & {a}_{j} & 0 & 0 & 1 & 1 & \\ & {b}_{j} & 0 & -1 & -1 & 0 & \\ & & & & & & \end{array}\end{eqnarray} \tag{ B8 }$$

such that random walk W_M (51a) pauses at first two 'time' moments and steps on two consecutive ones. The random walk V_M (51b) pauses at the first time moment, then makes two steps and pauses again.

Now, we again assume that M is divisible by 4. Then, the kernel in (52a) reduces to

$$\begin{eqnarray}\displaystyle \prod _{j=1}^{M}\cos ({a}_{j}x+{b}_{j}y) & = & {[\cos x\cos y\cos (x-y)]}^{\tfrac{M}{4}}\\ & = & \displaystyle \sum _{m=0}^{\tfrac{M}{4}}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{4}}{\tfrac{M}{4}-m}\right){(\cos x\cos y)}^{\tfrac{M}{2}-m}{(\sin x\sin y)}^{m},\end{eqnarray} \tag{ B9 }$$

such that the contributions dependent on x and on y appear in the factorized form. Inserting the latter expansion into equation (52a) and integrating term by term, we get

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\lambda }=\displaystyle \sum _{m=0}^{\tfrac{M}{4}}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{4}}{m}\right)\,{{ \mathcal U }}^{2}\left(\lambda ;\displaystyle \frac{M}{4}+m,\displaystyle \frac{M}{4}-m\right),\end{eqnarray} \tag{ B10 }$$

where ${ \mathcal U }$ is defined by

$$\begin{eqnarray}&&{ \mathcal U }(\lambda ;m,n)=\sqrt{\displaystyle \frac{M}{8\pi D{{\rm{\Delta }}}^{2}\lambda }}{\displaystyle \int }_{-\infty }^{+\infty }{\rm{d}}x\,{\cos }^{m}x{\sin }^{n}x\exp \left(-\displaystyle \frac{M}{8D{{\rm{\Delta }}}^{2}\lambda }{x}^{2}\right).\end{eqnarray} \tag{ B11 }$$

The function ${ \mathcal U }(\lambda ;m,n)$ is identically equal to zero for odd n. For even n, n = 2k, we write

$$\begin{eqnarray}&&{\cos }^{m}x{\sin }^{2k}x={\cos }^{m}x{(1-{\cos }^{2}x)}^{k}=\displaystyle \sum _{p=0}^{k}{(-1)}^{p}\left(\displaystyle \genfrac{}{}{0em}{}{k}{p}\right){\cos }^{2p+m}x.\end{eqnarray} \tag{ B12 }$$

Further on, for even m, m = 2q, the latter sum can be written down as a finite series of cosines

$$\begin{eqnarray}{\cos }^{2q}x{\sin }^{2k}x & = & {\cos }^{2q}x{(1-{\cos }^{2}x)}^{k}\\ & = & \displaystyle \sum _{p=0}^{k}\displaystyle \frac{{(-1)}^{p}}{{4}^{p+q}}\left(\displaystyle \genfrac{}{}{0em}{}{k}{p}\right)\displaystyle \sum _{l=-p-q}^{p+q}\left(\displaystyle \genfrac{}{}{0em}{}{2p+2q}{p+q+l}\right)\cos (2{lx}),\end{eqnarray} \tag{ B13 }$$

such that ${ \mathcal U }(\lambda ;2q,2k)$ is given explicitly by

$$\begin{eqnarray}&&{ \mathcal U }(\lambda ;2q,2k)=\displaystyle \sum _{p=0}^{k}\displaystyle \frac{{(-1)}^{p}}{{4}^{p+q}}\left(\displaystyle \genfrac{}{}{0em}{}{k}{p}\right)\displaystyle \sum _{l=-p-q}^{p+q}\left(\displaystyle \genfrac{}{}{0em}{}{2p+2q}{p+q+l}\right)\exp \left(-\displaystyle \frac{8D{{\rm{\Delta }}}^{2}{l}^{2}}{M}\lambda \right)\end{eqnarray} \tag{ B14 }$$

For odd m, m = 2q + 1, we have

$$\begin{eqnarray}{\cos }^{2q+1}x{\sin }^{2k}x & = & \displaystyle \frac{1}{2}\displaystyle \sum _{p=0}^{k}\displaystyle \frac{{(-1)}^{p}}{{4}^{p+q}}\left(\displaystyle \genfrac{}{}{0em}{}{k}{p}\right)\displaystyle \sum _{l=-p-q}^{p+q}\left(\displaystyle \genfrac{}{}{0em}{}{2p+2q}{p+q+l}\right)\\ & & \times (\cos ((2l+1)x)+\cos ((2l-1)x)),\end{eqnarray} \tag{ B15 }$$

which gives

$$\begin{eqnarray}{ \mathcal U }(\lambda ;2q+1,2k) & = & \displaystyle \frac{1}{2}\displaystyle \sum _{p=0}^{k}\displaystyle \frac{{(-1)}^{p}}{{4}^{p+q}}\left(\displaystyle \genfrac{}{}{0em}{}{k}{p}\right)\displaystyle \sum _{l=-p-q}^{p+q}\left(\displaystyle \genfrac{}{}{0em}{}{2p+2q}{p+q+l}\right)\\ & & \times \left(\exp \left(-\displaystyle \frac{8D{{\rm{\Delta }}}^{2}{(2l+1)}^{2}}{M}\lambda \right)\,+\,\exp \left(-\displaystyle \frac{8D{{\rm{\Delta }}}^{2}{(2l-1)}^{2}}{M}\lambda \right)\right).\end{eqnarray} \tag{ B16 }$$

Equations (B10), (B14) and (B16) define completely the moment-generating function of a single-trajectory periodogram in the case f = π/(2Δ). The distribution function follows straightforwardly by expanding ${{ \mathcal U }}^{2}$ in a series of exponentials and taking the inverse Laplace transform.

In figure B2 we plot our exact result for the moment-generating function in equation (B10) and the corresponding result for the distribution function of a single-trajectory periodogram, together with the conjectured expressions (56) and (57) with ${\mu }_{R}(f,T)$ and γ_R defined by equations (47a) and (47c). We observe a very convincing convergence of the exact results, upon an increase of M, to the conjectured forms.

Zoom In Zoom Out Reset image size

B.3. The case fΔ = 2π/3

In this case the coefficients a_j (42a) and b_j (42b) are periodic functions of j with period 3, and are given by

$$\begin{eqnarray}\begin{array}{cccccc} & & & & & \\ & j & 1 & 2 & 3 & \\ & {a}_{j} & 0 & 1/2 & 1 & \\ & {b}_{j} & 0 & -\sqrt{3}/2 & 0 & \\ & & & & & \end{array}\end{eqnarray} \tag{ B17 }$$

within their period.

We assume now that M is divisible by 3. In this case the kernel

$$\begin{eqnarray}\displaystyle \prod _{j=1}^{M}\cos ({a}_{j}x+{b}_{j}y) & = & {\left(\cos x\cos \left(\displaystyle \frac{x-\sqrt{3}y}{2}\right)\right)}^{\tfrac{M}{3}}={\cos }^{\tfrac{M}{3}}x\\ & & \times \displaystyle \sum _{p=0}^{\tfrac{M}{3}}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{3}}{p}\right){\left(\cos \displaystyle \frac{x}{2}\cos \displaystyle \frac{\sqrt{3}y}{2}\right)}^{p}{\left(\sin \displaystyle \frac{x}{2}\sin \displaystyle \frac{\sqrt{3}y}{2}\right)}^{\tfrac{M}{3}-p}.\end{eqnarray} \tag{ B18 }$$

Next, we use

$$\begin{eqnarray}&&{\cos }^{\tfrac{M}{3}}x={\left(2{\cos }^{2}\displaystyle \frac{x}{2}-1\right)}^{\tfrac{M}{3}}=\displaystyle \sum _{k=0}^{\tfrac{M}{3}}{(-1)}^{k}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{3}}{k}\right){\left(2{\cos }^{2}\displaystyle \frac{x}{2}\right)}^{\tfrac{M}{3}-k}\end{eqnarray} \tag{ B19 }$$

which gives, in combination with (B18), the following representation of the kernel

$$\begin{eqnarray}\displaystyle \prod _{j=1}^{M}\cos ({a}_{j}x+{b}_{j}y) & = & {2}^{\tfrac{M}{3}}\displaystyle \sum _{k=0}^{\tfrac{M}{3}}\displaystyle \frac{{(-1)}^{k}}{{2}^{k}}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{3}}{k}\right)\displaystyle \sum _{p=0}^{\tfrac{M}{3}}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{3}}{p}\right)\\ & & \times {\cos }^{p+\tfrac{2M}{3}-2k}\displaystyle \frac{x}{2}{\sin }^{\tfrac{M}{3}-p}\displaystyle \frac{x}{2}{\cos }^{p}\displaystyle \frac{\sqrt{3}y}{2}{\sin }^{\tfrac{M}{3}-p}\displaystyle \frac{\sqrt{3}y}{2}.\end{eqnarray} \tag{ B20 }$$

Inserting the latter expansion into equation (52a) and integrating term by term, we obtain the following exact result for the moment-generating function of a single-trajectory periodogram

$$\begin{eqnarray}{{\rm{\Phi }}}_{\lambda } & = & {2}^{\tfrac{M}{3}}\displaystyle \sum _{k=0}^{\tfrac{M}{3}}\displaystyle \frac{{(-1)}^{k}}{{2}^{k}}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{3}}{k}\right)\displaystyle \sum _{p=0}^{\tfrac{M}{3}}\left(\displaystyle \genfrac{}{}{0em}{}{\tfrac{M}{3}}{p}\right)\\ & & \times { \mathcal U }\left(\displaystyle \frac{\lambda }{4};p+\displaystyle \frac{2M}{3}-2k,\displaystyle \frac{M}{3}-p\right)\,{ \mathcal U }\left(\displaystyle \frac{3\lambda }{4};p,\displaystyle \frac{M}{3}-p\right),\end{eqnarray} \tag{ B21 }$$

where ${ \mathcal U }$ is defined in equations (B14) and (B16). The moment-generating function is a finite sum of exponentials and the corresponding distribution is a sum of delta-functions with coefficients which can be readily deduced from equation (B21). In figure B3 we demonstrate that Φ_λ in equation (B21), and also the cumulative distribution function of the distribution deduced from equation (B21), do indeed converge to the conjectured forms (56) and (57) with ${\mu }_{R}(2\pi /(3{\rm{\Delta }}),T)$ and ${\sigma }_{R}^{2}(2\pi /(3{\rm{\Delta }}),T)$, equations (48a) and (48b).

Zoom In Zoom Out Reset image size

B.4. The case fΔ = 2π