<figure><img src="image_1.png" /></figure>

HAL
open science

# $1/f^{3/2}$ noise in heavy traffic M/M/1 queue

Thomas Vanderbruggen

► To cite this version:

Thomas Vanderbruggen. $1/f^{3/2}$ noise in heavy traffic M/M/1 queue. 2023. hal-03973799v2

HAL Id: hal-03973799

https://hal.science/hal-03973799v2

Preprint submitted on 22 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
$1/f^{3/2}$ NOISE IN HEAVY TRAFFIC M/M/1 QUEUE

T. VANDERBRUGGEN

**ABSTRACT.** We study the power spectral density of continuous time Markov chains and explicit its relationship with the eigenstructure of the infinitesimal generator. This result helps us understand the dynamics of the number of customers for a M/M/1 queuing process in the heavy traffic regime. Closed-form relations for the power law scalings associated to the eigenspectrum of the M/M/1 queue generator are obtained, providing a detailed description of the power spectral density structure, which is shown to exhibit a $1/f^{3/2}$ noise. We confirm this result by numerical simulation. We also show that a continuous time random walk on a ring exhibits similar behavior with the same scaling exponent. It is remarkable that a complex behavior such as $1/f$ noise can emerge from the M/M/1 queue, which is the "simplest" queuing model.

# 1. INTRODUCTION

Flicker noise also named $1/f$ noise [28, 40] describes a stochastic process whose power spectral density (PSD) exhibits a power law behavior $1/f^\alpha$ over a certain range of frequencies $f$, usually dominating at low frequency. The cases $\alpha = 0$ and $\alpha = 2$ correspond to white noise and Brownian motion noise, respectively. Both are well understood and are simply related to each other by derivation / integration - The Brownian motion being the integral of the white noise. Hence a $1/f^\alpha$ noise is usually defined for $\alpha$ in the range $0 < \alpha < 2$, excluding the trivial cases $\alpha = 0$ and $\alpha = 2$. In other words a $1/f^\alpha$ noise cannot be obtained directly from simple transformation on the white noise, that is the reason why its mathematical description is challenging.

On the experimental side the flicker noise appears in a wide range of systems. For example: heart rate [1], human cognition [11], music [38], internet traffic [6, 34], electronic devices [35], magnetic field fluctuations of the Sun [13] and Saturn [14]. However, it is still an open question whether the widespread existence of this behavior is explained by a common underlying mathematical mechanism.

Several mechanisms have been proposed to generate a $1/f^\alpha$ noise, and a detailed presentation can be found in [40]. The two main mathematical frameworks are the superposition of Poisson processes with widespread relaxation rates [39, 44, 22, 8, 9], and non-linear stochastic differential equations [32, 33]. The superposition of relaxation rates, whose resulting PSD is a sum of Lorentzians, is an interesting candidate, especially because reversible homogeneous continuous time Markov processes exhibit such a spectrum, as shown in [22, 8]. Here, we redemonstrate this result and show that the cut-off frequencies are the eigenvalues of the infinitesimal generator and the amplitude of each Lorentzian is given by the coupling of the state vector to the eigenvectors. In particular, we show that if the eigenvalues and the eigenvectors follow appropriate power laws a $1/f^\alpha$ noise scaling can be achieved. However, to the extent of our knowledge, mostly *ad hoc* Markov chain with tailored eigenspectrum distributions have been introduced to show that such models can generate $1/f$ noise. Here we show that $1/f$ noise appears in a "classical" Markov model, namely the M/M/1 queue.

The M/M/1 queue is a birth-death process with constant transition rates, making it the simplest birth-death processes (beyond pure birth and pure death processes). A major advantage of the M/M/1 queue simplicity is that many closed-form expressions are known for various queue metrics (see for example [15]). We apply the previously introduced result on the PSD of reversible Markov chains to the infinitesimal generator of a M/M/1 queue in the heavy traffic regime. Scaling laws for the eigenvalues and eigenvectors of the generator are obtained from its diagonalization and used to show that the PSD exhibits a $1/f^{3/2}$ noise, which is confirmed by numerical simulation. Finally, we consider a continuous time random walk on a ring and show that its generator eigenstructure is similar to the M/M/1 queue generator, and therefore it also exhibits a $1/f^{3/2}$ noise as was already noticed [8, 43].

<footnote>Licence CC-BY 4.0
Date: October 22, 2023.</footnote>
## 2. DEFINITIONS AND NOTATIONS

**2.1. Continuous time Markov chains.** We consider homogeneous continuous time Markov chains (HCTMC) on a finite countable state space of size $n$. A HCTMC $X = \{X_t, t \ge 0\}$ is characterized by its transition matrix $P(\tau)$ with elements $P_{ij}(\tau) = P(X_\tau = j|X_0 = i)$. It is a stochastic matrix, because $\sum_{i=1}^n P_{ij}(\tau) = 1$, which obeys to the Kolmogorov forward equation:

$$ (1) \qquad \frac{d}{d\tau} \mathbf{P}(\tau) + \mathbf{G} \mathbf{P}(\tau) = \mathbf{0}, $$

where $\mathbf{G}$ is called the infinitesimal generator of the HCTMC. Since a continuous time Markov chain over a finite state space is regular, there is exactly one solution to this differential equation. Given the initial condition $\mathbf{P}(0) = \mathbf{1}$ must be satisfied, the solution of this differential equation is

$$ (2) \qquad \mathbf{P}(\tau) = e^{-\mathbf{G}\tau}. $$

The generator satisfies $\sum_i \mathbf{G}_{ij} = 0$ and all the off-diagonal elements must be negative.

The long term behavior of the Markov process is characterized by the stationary probability vector $\pi = (\pi_i)_{1 \le i \le n}$, where $\pi_i$ is the probability to be in the state $i$ of the chain after an infinite time. We further assume the chain is irreducible so that the stationary distribution is unique, because the state space is finite. From Eq. (1), it satisfies $\pi \mathbf{G} = \mathbf{0}$ with $\|\pi\|_1 = 1$.

Finally, we consider a reversible chain, that is a chain which satisfies the detailed balance equations $\pi_i \mathbf{G}_{ij} = \pi_j \mathbf{G}_{ji}$. A reversible Markov process has the desirable property to be stationary, and therefore the power spectral density of this process exists. Because the matrix $\mathbf{D}^{1/2} \mathbf{G} \mathbf{D}^{-1/2}$, where $\mathbf{D} = \text{diag}(\pi)$, is symmetric and similar to $\mathbf{G}$, the generator of a reversible chain is diagonalizable with real positive eigenvalues [18]. Note that birth-death processes, and in particular the M/M/1 queue, are reversible (see for example [15] 8.4.3).

**2.2. Scalar product.** In the context of Markov chains the scalar product whose metric is the stationary probability occurs naturally. To each state $i$ of the chain we associate a real value $x_i$, and we note $\mathbf{x} = (x_i)_{1 \le i \le n} \in \mathbb{R}^n$. For a HCTMC $X$ with stationary probability $\pi$, the space $\mathbb{R}^n$ is equipped with the scalar product

$$ (3) \qquad \langle \mathbf{u}, \mathbf{v} \rangle_{\pi} \triangleq \sum_{i=1}^{n} \pi_i u_i v_i. $$

The all-ones vector is noted $\mathbf{1}$, it verifies $\langle \mathbf{1}, \mathbf{1} \rangle_{\pi} = \|\pi\|_1 = 1$. This scalar product relates to statistical quantities such as the average or second order moment:

$$ (4) \qquad \langle \mathbf{x}, \mathbf{1} \rangle_{\pi} = \sum_{i=1}^{n} \pi_i x_i = \langle \mathbf{x} \rangle, $$

$$ (5) \qquad \langle \mathbf{x}, \mathbf{x} \rangle_{\pi} = \sum_{i=1}^{n} \pi_i x_i^2 = \langle x^2 \rangle. $$

## 3. POWER SPECTRAL DENSITY OF A MARKOV PROCESS

We re-establish a known result that the PSD of a reversible Markov process is a superposition of Lorentzians [22, 8]. The derivation explicits the relationship with the fundamental matrix of the Markov chain as an intermediate result. We first derive the expression for the autocorrelation function of a stationary Markov process, and, combining this result with the Green-Kubo relation, we show the relation between the diffusion coefficient of the process and the fundamental matrix of the Markov chain. Interestingly, the Green-Kubo relation is somehow a particular case of the Wiener-Kinchine theorem (the diffusion coefficient is the PSD at zero frequency), and we show that the PSD of the Markov chain relates to a generalized version of the fundamental matrix. Then we explicit the relationship between the generalized fundamental matrix and the infinitesimal generator, and we show that the PSD of a reversible process is a sum of Lorentzians. Finally, we establish a criterion relating power law scalings of the generator eigenstructure and $1/f$ noise scaling.
**3.1. Autocorrelation function.** Here we consider the centered process $\delta X_t \triangleq X_t - \langle X_t \rangle$, in other word we remove the DC component from the power spectrum. The autocorrelation function is defined as

$$ (6) \qquad C_X(\tau) \triangleq \lim_{t \to \infty} \langle \delta X_{t+\tau} \delta X_t \rangle. $$

In App. A we prove that the autocorrelation function of a HCTMC is a quadratic form:

**Proposition 1.** Let $X$ be a stationary and irreducible HCTMC and let $\mathbf{x}$ be a real-valued vector over the (finite) state space of the chain, then the autocorrelation of $X$ is

$$ (7) \qquad C_X(\tau) = \langle \mathbf{x}, [\mathbf{P}(\tau) - \mathbf{P}(\infty)] \mathbf{x} \rangle_{\pi}. $$

Where we defined the limiting transition matrix $\mathbf{P}(\infty)$, whose lines are the stationary probability vector $\mathbf{P}_{ij}(\infty) \triangleq \lim_{\tau \to \infty} \mathbf{P}_{ij}(\tau) = \pi_j$.

Note that this result and the next ones assume a finite state space with implies that the stationary population is unique, but the results would hold if an unique distribution exists, whether the state space is finite or not.

**3.2. Fundamental matrix and Green-Kubo relation.** The fundamental matrix, also called deviation matrix, of a HCTMC is defined as

$$ (8) \qquad \mathbf{Z} \triangleq \int_{0}^{\infty} [\mathbf{P}(\tau) - \mathbf{P}(\infty)] d\tau. $$

It can be related to many properties of the Markov chain, such as variance [41], first-passage time analysis [45], speed of convergence to the stationary distribution [5]. Interestingly, we see here that, by the Green-Kubo relation, it is also related to the diffusion coefficient of the process, and more generally to some transport coefficient of the process.

The diffusion coefficient of the process $X$ is $D_X = D_{\partial_t Y} \triangleq \lim_{\tau \to \infty} \sigma_Y^2(t)/2t$, where $\sigma_Y^2(t)$ is the variance of the process $Y$ which is the integral of the process $X$ ($X = \partial_t Y$) [10]. The diffusion coefficient is related to the autocorrelation function by the Green-Kubo relation [7, 10]:

$$ (9) \qquad D_X = \int_0^\infty C_X(\tau) d\tau. $$

Using the result Prop. 1 for the autocorrelation function, it follows

$$ (10) \qquad D_X = \left\langle \mathbf{x}, \left\{ \int_0^\infty [\mathbf{P}(\tau) - \mathbf{P}(\infty)] d\tau \right\} \mathbf{x} \right\rangle_{\pi}, $$

and we observe that the Green-Kubo relation for the diffusion coefficient of a stationary Markov process is a quadratic form of the fundamental matrix:

**Proposition 2 (Green-Kubo relation).** Let $X$ be a stationary and irreducible HCTMC and $\mathbf{x}$ be a real-valued vector indexed over the (finite) state space of the chain, then the diffusion coefficient of $X$ is

$$ (11) \qquad D_X = \langle \mathbf{x}, \mathbf{Z}\mathbf{x} \rangle_{\pi}. $$

**3.3. Power spectral density.** The Wiener-Kinchine theorem states that the power spectral density is the Fourier transform of the autocorrelation function. More precisely, if $X$ is a wide-sense stationary process then the power spectral density of $X$ is

$$ (12) \qquad S_X(\omega) = \int_0^\infty C_X(\tau) \cos(\omega\tau) d\tau. $$

We note that this is related to the Green-Kubo relation Eq. (9) because $D_X = S_X(0)$, and the power spectral density also describes out-of-equilibrium fluctuations around the stationary state of the system.

Because a stationary process is wide-sense stationary, applying the Wiener-Kinchine theorem to the result Prop. 1 for the autocorrelation function yield:

$$ (13) \qquad S_X(\omega) = \left\langle \mathbf{x}, \left\{ \int_0^\infty [\mathbf{P}(\tau) - \mathbf{P}(\infty)] \cos(\omega\tau) d\tau \right\} \mathbf{x} \right\rangle_{\pi}. $$
Now we define a generalized version of the fundamental matrix Eq. (8):

$$ (14) \qquad \mathbf{Z}(\omega) \triangleq \int_{0}^{\infty} [\mathbf{P}(\tau) - \mathbf{P}(\infty)] \cos(\omega\tau) \, d\tau, $$

where the fundamental matrix is retrieved at the zero frequency: $\mathbf{Z} = \mathbf{Z}(0)$. Hence the Wiener-Kinchine theorem for a stationary Markov chain states that the PSD is a quadratic form with matrix $\mathbf{Z}(\omega)$:

**Proposition 3** (Wiener-Kinchine). Let $X$ be a stationary and irreducible HCTMC and $\mathbf{x}$ be a real-valued vector indexed over the (finite) state space of the chain, then the power spectral density of $X$ is

$$ (15) \qquad S_X(\omega) = \langle \mathbf{x}, \mathbf{Z}(\omega)\mathbf{x} \rangle_{\pi}. $$

The following property, proven in App. B not assuming $\mathbf{G}$ to be diagonalizable, relates the fundamental matrix to the infinitesimal generator:

**Proposition 4.** The generalized fundamental matrix of an irreducible HCTMC over a finite state space with generator $\mathbf{G}$ is

$$ (16) \qquad \mathbf{Z}(\omega) = (\mathbf{G} + \omega^2 \mathbf{G}^{\sharp})^{\sharp}, $$

where $\mathbf{A}^{\sharp}$ is the group inverse of the matrix $\mathbf{A}$.

In particular, $\mathbf{Z} = \mathbf{Z}(0) = \mathbf{G}^{\sharp}$ and we recover the known fact that the fundamental matrix is the group inverse of the generator (see [5] for example). The Green-Kubo relation Eq. (11) is thus $D_X = \langle \mathbf{x}, \mathbf{G}^{\sharp}\mathbf{x} \rangle_{\pi}$, as was observed in [36].

3.4. **Sum of Lorentzians.** In the case of a reversible chain, which is diagonalizable with real positive eigenvalues, we can combine the results of Prop. 3 and Prop. 4 to obtain an explicit expression relating the power spectral density to eigenspectrum structure of the generator (see details in App. C):

**Theorem 5.** Let $X$ be an irreducible and reversible HCTMC over a finite state space of size $n$. If the non-zero eigenvalues of the infinitesimal generator of $X$ are $\{\omega_k\}_{1\le k\le n-1}$, and $\mathbf{\Pi}_k$ is the spectral projector associated to the eigenvalue $\omega_k$, then the power spectral density of $X$ is

$$ (17) \qquad S_X(\omega) = \sum_{k=1}^{n-1} \gamma_k^2 \frac{\omega_k}{\omega_k^2 + \omega^2}, $$

where $\gamma_k^2 \triangleq \langle \mathbf{x}, \mathbf{\Pi}_k \mathbf{x} \rangle_{\pi}$.

We observe that the power spectral density of a Markov chain is a sum of Lorentzian. Similar results are obtained for reversible Markov chain [22, 8, 9], nonlinear stochastic differential equations [32] or the spectral density of the laser field near threshold [27, 18.7.3]. Note that since the process is reversible every non-zero eigenvalue of $\mathbf{G}$ is real and strictly positive, so that $S_X(\omega) \ge 0$ as expected.

The coefficient $\gamma_k^2$ corresponds to the coupling strength of the $k$th eigenvector to the state vector $\mathbf{x}$. In App. F, we show that for symmetric infinitesimal generators $\gamma_k$ is related to the graph Fourier transform of the signal $\mathbf{x}$ over the vertices of the (undirected) graph associated to the Markov chain.

**Remark 1** (Normalization of $\gamma_k$ coefficients). Since the spectral projectors satisfies the resolution of unity $\sum_{k=1}^{n-1} \mathbf{\Pi}_k = \mathbf{1}$, the coupling strength has the normalization:

$$ (18) \qquad \sum_{k=1}^{n-1} \gamma_k^2 = \langle \mathbf{x}, \mathbf{x} \rangle_{\pi} = \langle x^2 \rangle. $$

**Remark 2** (Normalization of $S_X(\omega)$). The overall energy of the spectrum is, using Eq. (18):

$$ (19) \qquad \int_{0}^{\infty} S_X(\omega) \, d\omega = \sum_{k=1}^{n-1} \gamma_k^2 \int_{0}^{\infty} \frac{\omega_k}{\omega_k^2 + \omega^2} \, d\omega = \frac{\pi}{2} \sum_{k=1}^{n-1} \gamma_k^2 = \frac{\pi}{2} \langle x^2 \rangle. $$

Therefore if we desire to satisfy $\int_{0}^{\infty} S_X(\omega) \, d\omega = \langle x^2 \rangle$, we must include a normalization factor $2/\pi$ to the result Eq. (17).
**3.5. Lorentzian power spectrum.** In the specific case where the eigenspectrum of **G** is degenerate ($\omega_k = \bar{\omega}$ for all $k > 0$), which occurs in particular for a two states chain, then the power spectrum is Lorentzian:

$$ (20) \qquad S_X(\omega) = \frac{\langle x^2 \rangle}{\bar{\omega}} \frac{1}{1 + (\omega/\bar{\omega})^2} = \frac{D_X}{1 + (\omega/\bar{\omega})^2}. $$

*Remark 3.* Another, more contrived, way to obtain a Lorentzian spectrum is for **x** to be an eigenvector of **G**.

As an example of a process with a highly degenerate spectrum we consider the continuous time random walk on star graph. It is a process where one transition from a central state to a peripherical one with rate $\lambda$, and back to the central state from the periphery with rate $\mu$. The transition graph is presented in Fig. 1.

<figure><img src="image_3.png" /><figcaption>FIGURE 1. Transitions for the random walk on a star graph.</figcaption></figure>

The infinitesimal generator of this process is an arrowhead matrix:

$$ (21) \qquad \mathbf{G} = \begin{pmatrix} (n-1)\lambda & -\lambda & -\lambda & \dots \\ -\mu & \mu & & \\ -\mu & & \mu & \\ \vdots & & & \vdots \end{pmatrix}. $$

We further assume $\lambda = \mu = 1$ and the generator is the Laplacian matrix of an undirected star graph. The eigenvalues are $\omega_0 = 0$, $\omega_1 = \omega_2 = \dots = \omega_{n-2} = 1$ and $\omega_{n-1} = n$ (see for example [3] 15.3.3 (ii)). Therefore, excluding the zero eigenvalue, the generator only has 2 eigenvalues 1 and n so the spectrum is strongly degenerate and the power spectral density is nearly Lorentzian, as we can observe on the simulation Fig. 2.

*Remark 4.* In the special case where $n = 1$ one obtain the telegraph process which gives rise to random telegraph noise (also known as burst noise or popcorn noise), whose power spectrum is well known to be Lorentzian [26].

**3.6. 1/f noise criterion.** When the eigenspectrum of **G** is highly non-degenerate, the power spectral density can present more complex shapes beyond the Lorentzian spectrum. In fact such a spectrum, resulting from the superposition of relaxation processes, has often been proposed as an explanation for 1/f noises [28]. It often assumes an uniform distribution of the eigenfrequencies over a given range. More generally, as already observed by [9], the power spectrum is determined by the complete eigenstructure ($\omega_k, \gamma_k$) of the generator. In particular, if $\omega_k$ and $\gamma_k$ obey appropriate power laws then the power spectral density exhibits a 1/f like noise, more precisely the criterion is:

**Criterion 6.** Assume that the eigenvalues of **G** scales as $\omega_k \sim k^\alpha$ with $\alpha > 0$, and that the coupling strength of the eigenvectors to the state vector **x** scales as $\gamma_k \sim k^\beta$.

If $\alpha > |2\beta + 1|$, then the power spectral density scales as $S_X(\omega) \sim \omega^\zeta$, where

$$ (22) \qquad \zeta = \frac{2\beta - \alpha + 1}{\alpha}. $$
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Simulation of the continuous time random walk on a star graph of size $n = 1000$. (a) Example of simulation trajectory (state on the ring versus time). (b) Periodogram (averaged over 20 realizations), the green line is a fit with a Lorentzian function.</figcaption></figure>

It satisfies $-2 < \zeta < 0$. An exact $1/f$ noise ($\zeta = -1$) is obtained for $\beta = -1/2$, for any $\alpha > 0$.

In practice the scaling law may only be valid over a finite range of eigenfrequencies and to achieve a $1/f^\alpha$ noise spectrum the eigenspectrum of $G$ must be highly degenerate with eigenvalues spanning over at least one frequency decade. As we will see now, this property appears for queues (birth-death processes) in the heavy traffic regime.

## 4. APPLICATION TO THE M/M/1 QUEUE IN THE HEAVY TRAFFIC REGIME

A M/M/1 queue is a birth-death process with constant rates ($\mu = \mu_n$, $\lambda = \lambda_n$ for all $n$). The associated graph is presented in Fig. 3. It models a queue with a single server where arrival time and server processing time obey an exponential distribution with time constants independent of the number of customers in the queue. This is the "simplest" queue model and many closed-form results are known.

<figure><img src="image_5.png" /><figcaption>FIGURE 3. Transition graph of the M/M/1 queue.</figcaption></figure>

The infinitesimal generator of this process is a tridiagonal matrix:

$$ (23) \qquad \mathbf{G} = \begin{pmatrix} \lambda & -\lambda & & & \\ -\mu & \lambda+\mu & -\lambda & & \\ & -\mu & \lambda+\mu & -\lambda & \\ & & \ddots & \ddots & \\ & & & \lambda & \end{pmatrix}, $$

The server utilization is defined as $\rho \triangleq \lambda/\mu$. The process is stable when $\rho < 1$, that is when the long term number of customers in the queue is finite. The average number of customers in the queue is $m = \rho/(1-\rho)$ and the variance is $\sigma^2 = \rho/(1-\rho)^2$. The Markov chain has a stationary population which satisfies $\pi_k = (1-\rho) \rho^k$.
4.1. **Tridiagonal Toeplitz matrix.** The tridiagonal Toeplitz matrix is a $n \times n$ matrix with constant upper diagonal, lower diagonal and diagonal:

$$ (24) \qquad T_n(a, b, c) \triangleq \begin{pmatrix} b & c & & & \\ a & b & c & & \\ & a & \ddots & \ddots & \\ & & \ddots & b & c \\ & & & a & b \end{pmatrix}, $$

In the limit $n \to \infty$, the generator Eq. (23) can be approximated by $T_n(a, b, c)$, where $a = -\mu$, $b = \mu + \lambda$, and $c = -\lambda$. The only difference with the M/M/1 generator being the first row does not satisfies the condition $b + c = 0$.

The eigenvalues and eigenvectors of $T_n(a, b, c)$ have simple closed-form expressions. The eigenvalues (indexed in increasing order) in the case $a, c < 0$ are (see for example [30]):

$$ (25) \qquad \omega_k = b - 2\sqrt{ac} \cos \left( \frac{k\pi}{n+1} \right). $$

The right eigenvectors are $v_k = (v_{k,i})_{1 \le i \le n}$ where

$$ (26) \qquad v_{k,i} = \left(\frac{a}{c}\right)^{i/2} \sin\left(\frac{ki\pi}{n+1}\right), $$

and the left eigenvectors $w_k = (w_{k,i})_{1 \le i \le n}$ are:

$$ (27) \qquad w_{k,i} = \left(\frac{c}{a}\right)^{i/2} \sin\left(\frac{ki\pi}{n+1}\right). $$

When the matrix is symmetric ($c = a$), the eigenvectors are the discrete Fourier transform basis, this is because the infinitesimal generator is the Laplacian matrix of the graph associated to the chain. In fact in App. F we explicit the relation between the power spectral density and the graph Fourier transform. In the case of a birth-death process, the associated graph is a discrete semi-open unidimensional string whose eigenmodes (the eigenvectors of the graph Laplacian matrix) are the discrete Fourier transform basis.

4.2. **Light traffic regime.** In the light traffic regime, the server utilization is low $\rho \ll 1$, and there are few customers in the queue. Choosing $\mu = 1$ and $\lambda = \varepsilon$, the server utilization is $\rho = \varepsilon$. The average number of customers is $m \sim \varepsilon$ and the variance is $\sigma^2 \sim \varepsilon$, in other words $\sigma \sim \sqrt{m}$, and the underlying distribution behaves as a Poisson distribution. In the limit $\varepsilon \to 0$, the eigenvalues satisfy

$$ (28) \qquad \omega_k = 1 + \varepsilon - 2\sqrt{\varepsilon} \cos \left( \frac{k\pi}{n+1} \right) \underset{\varepsilon \to 0}{\longrightarrow} 1, $$

so that all the eigenfrequencies are degenerate and the power spectral density is Lorentzian.

4.3. **Scaling laws in the heavy traffic regime.** In the heavy traffic regime the server utilization is high ($\rho \sim 1$) and the number of customers in the queue is large. This regime is known to converge to a reflected Brownian motion [23, 21], that is a Brownian motion over the half-plane $\mathbb{R}^+$. Here, we set $\lambda = 1$ and $\mu = \lambda + \varepsilon = 1 + \varepsilon$ where $\varepsilon \ll 1$. Therefore the utilization is $\rho \sim 1 - \varepsilon$, the average number of customers is $m \sim 1/\varepsilon$ and the variance is $\sigma^2 \sim 1/\varepsilon^2$. So in this regime the average number of customers equals the standard deviation ($\sigma \approx m$), showing that the underlying distribution has much larger fluctuations than a Poisson process.

From Eq. (25), the eigenvalues of the generator satisfies:

$$ (29) \qquad \omega_k = 2 + \varepsilon - 2\sqrt{1+\varepsilon} \cos \left( \frac{k\pi}{n+1} \right) \underset{\varepsilon \to 0}{\sim} 2 \left[ 1 - \cos \left( \frac{k\pi}{n} \right) \right] \underset{k \ll n}{\sim} \left( \frac{k\pi}{n} \right)^2. $$

In other words the eigenfrequencies follow the power law $\omega_k \sim k^2$ ($\alpha = 2$).

In Fig. 4 (a), we compare this approximation with the exact analytical result for the Toeplitz matrix, and with the numerical eigenvalues for the M/M/1 generator. We observe that for large $n$, the eigenfrequencies are highly non-degenerate and are densely covering many order of magnitudes. In appendix E, we see that this behavior is generalized among birth-death processes because the
characteristic polynomial of a tridiagonal matrix is an orthogonal polynomial whose roots are
densely spanning an interval.

<figure><img src="image_6.png" /><figcaption>FIGURE 4. Power law scalings for the M/M/1 infinitesimal generator. Evolution of (a) the eigenvalues $\omega_k$ and (b) the coefficients $|\gamma_k|$ for the generator and the tridiagonal Toeplitz matrix versus $k$ for $n = 10^3$, $\varepsilon = 10^{-4}$. The blue dots are the result of the numerical eigendecomposition of the M/M/1 queue generator, the green dots are the analytical result for the Toeplitz matrix and the solid line are the power law approximations $\omega_k \sim (k\pi/n)^2$ and $|\gamma_k| \sim \sqrt{\varepsilon n^2/(\pi k)}$.</figcaption></figure>

Since we analyze the number of customers in the queue versus time, the weight associated
to state $i$ is simply $x_i = i$. The coupling coefficient $\gamma_k$ satisfies (using the spectral projector
$\Pi_k = \mathbf{w}_k^\top \mathbf{v}_k$):

$$
(30) \qquad \gamma_k^2 = \langle \mathbf{x}, \mathbf{w}_k^\top \mathbf{v}_k \mathbf{x} \rangle_{\pi}
$$

$$
(31) \qquad = (1-\rho) \sum_{i,j} ij\rho^{(i+j)/2} \sin\left(\frac{ki\pi}{n+1}\right) \sin\left(\frac{kj\pi}{n+1}\right)
$$

$$
(32) \qquad = (1 - \rho) \left[ \sum_{i} i \rho^{i/2} \sin \left( \frac{ki\pi}{n+1} \right) \right]^2 .
$$

In the limit $\varepsilon \to 0$, $1 - \rho = \varepsilon + o(\varepsilon^2)$ and $\rho^{i/2} = 1 - i\varepsilon/2 + o(\varepsilon^2)$, hence

$$
(33) \qquad \gamma_k^2 = \varepsilon \left[ \sum_{i=1}^{n} i \sin \left( \frac{ki\pi}{n+1} \right) \right]^2 + o(\varepsilon^2).
$$

We can evaluate the sum:

$$
(34) \quad \sum_{i=1}^{n} i \sin \left( \frac{ki\pi}{n+1} \right) = - \frac{(n+1)}{2} \frac{\cos \left( \frac{\pi k (2n+1)}{2n+2} \right)}{\sin \left( \frac{\pi k}{2n+2} \right)} \underset{n \to \infty}{\sim} - \frac{n \cos(\pi k)}{2 \sin(\frac{\pi k}{2n})} \underset{k \ll n}{\sim} (-1)^{k+1} \frac{n^2}{\pi k}.
$$

Therefore $|\gamma_k| \sim \sqrt{\varepsilon n^2/(\pi k)}$ and $\gamma_k$ scales as $k^{-1}$ ($\beta = -1$). In Fig. 4 (b), we compare this approximation with the exact analytical result.

The scaling exponents $(\alpha, \beta) = (2, -1)$ satisfy $\alpha = 2 > |2\beta + 1| = 1$, and from Criterion 6 the power spectral density scales as $S_X(\omega) \sim \omega^\zeta$, where $\zeta = (2\beta - \alpha + 1)/\alpha = -3/2$. In other words,
the power spectral density of the number of customers in a M/M/1 queue in the heavy traffic regime exhibits a $1/f^{3/2}$ noise.

**4.4. Heavy traffic simulations.** The M/M/1 queue is simulated using the Gillepsie algorithm [12], where each iteration compute the remaining time in a given state and the jump to the next step. Since the remaining time in a given state is random, the Gillepsie algorithm does not produce uniformly sampled data, so the result is interpolated to obtain the uniform sampling required by the Fast Fourier Transform (FFT) algorithm. The power spectral density is estimated from the periodogram, obtained by computing the modulus square of the FFT of the data. Periodograms of many simulation trajectories are averaged to improve the power spectral density estimation.

<figure><img src="image_7.png" /><figcaption>FIGURE 5. Simulation of the M/M/1 queue in the heavy traffic regime. Example of simulation trajectories (number of customers in the queue versus time) for (a) $\epsilon = 10^{-3}$ and (c) $\epsilon = 10^{-4}$. Periodograms for (b) $\epsilon = 10^{-3}$ (averaged over 69 realizations) and (d) $\epsilon = 10^{-4}$ (averaged over 34 realizations).</figcaption></figure>

Simulation results are presented Fig. 5 for $\epsilon = 10^{-3}$ and $\epsilon = 10^{-4}$. We check that the mean value satisfies $m \sim 1/\epsilon$: for $\epsilon = 10^{-3}$, $m = 987$ (expected $10^3$), and for $\epsilon = 10^{-4}$, $m = 10357$ (expected $10^4$). We observe the periodogram scales as $1/f^{3/2}$ over at least a frequency decade. Note also the $1/f^2$ scaling at higher frequencies, showing the spectrum is consistent with a sum of Lorentzians.

## 5. CONTINUOUS TIME RANDOM WALK ON A RING

To further understand the above result, that is a continuous time random walk on a semi-open string, we fold the string onto itself to form a ring and consider the resulting continuous time random walk on a ring. This process could for example model the phase diffusion of an oscillator.

Interestingly, a $1/f^{3/2}$ noise has also been observed for a random walk on a ring [8, 43]. Here we confirm this result and show that the PSD of a continuous time random walk on a ring obeys a $1/f^{3/2}$ scaling. Indeed, we observe that the eigenstructure of the circular random walk generator
is similar with the one of the M/M/1 queue generator, more specifically the scaling law coefficients $(\alpha, \beta)$ are the same.

<figure><img src="image_8.png" /><figcaption>FIGURE 6. Transition graph of the continuous time random walk on a ring.</figcaption></figure>

5.1. **Generator, eigenvalues and eigenvectors.** We consider a ring of $n$ states with constant transition rates clockwise $\lambda$ and counterclockwise $\mu$ (transition graph is presented in Fig. 6). The infinitesimal generator of this process is a $n \times n$ matrix:

$$ (35) \qquad \mathbf{G} = \begin{pmatrix} \lambda + \mu & -\lambda & & & -\mu \\ -\mu & \lambda + \mu & -\lambda & & \\ & -\mu & \lambda + \mu & \ddots & \\ & & \ddots & \ddots & -\lambda \\ & & & -\mu & \lambda + \mu \end{pmatrix}. $$

It is a right circulant matrix, which is strongly related to the discrete Fourier transform. In what follows, we define $w_n \triangleq \exp(i2\pi/n)$ so that $w_n^k$ is the $k$th root of unity. The eigenvalues of the generator are:

$$ (36) \qquad \omega_k = \lambda + \mu - \lambda w_n^k - \mu w_n^{(n-1)k}. $$

We further assume that clockwise and counterclockwise circulations are identical, that is $\lambda = \mu \triangleq 1$. In that case the generator is the Laplacian matrix of an undirected cycle graph $C_n$. As a consequence the generator is a symmetric matrix and the Markov process is reversible. The eigenvalues satisfy:

$$ (37) \qquad \omega_k = 2 - 2 \cos \left( \frac{2\pi k}{n} \right) \underset{k \le n}{\sim} 4 \left( \frac{\pi k}{n} \right)^2. $$

The eigenvalues scaling exponent is therefore $\alpha = 2$. This result is very similar with the M/M/1 case Eq. (29), up to a factor 2 in the cosine argument.

The eigenvectors $\mathbf{v}_k = (v_{k,q})_{1 \le q \le n}$ are the same for any circulant matrix, in particular it does not depend on $\lambda$ and $\mu$:

$$ (38) \qquad v_{k,q} = \frac{1}{\sqrt{n}} w_n^{(q-1)k}. $$

Because the generator is symmetric, the right and left eigenvectors are identical. The spectral projector is $\mathbf{\Pi}_k = \mathbf{v}_k^\dagger \mathbf{v}_k$ where $\mathbf{v}_k^\dagger$ is the hermitian conjugate of $\mathbf{v}_k$, it satisfies:

$$ (39) \qquad (\mathbf{\Pi}_k)_{q,r} = \frac{1}{n} w_n^{(1-q)k} w_n^{(r-1)k} = \frac{1}{n} w_n^{(r-q)k}. $$
5.2. 1/f³/² **noise**. The coupling coefficient $\gamma_k^2$ for $x_q = q$ is:

$$ (40) \qquad \gamma_k^2 = \sum_{q,r} \pi_q q r (\Pi_k)_{q,r} $$

$$ (41) \qquad = \frac{1}{n^2} \sum_{q,r} q r w_n^{(r-q)k} $$

$$ (42) \qquad = \frac{1}{n^2} \left( \sum_{q=0}^{n-1} q w_n^{-qk} \right) \left( \sum_{r=0}^{n-1} r w_n^{rk} \right) $$

$$ (43) \qquad = \frac{1}{n^2} \frac{n}{2ie^{i\pi k/n} \sin(\pi k/n)} \frac{ne^{i\pi k/n}}{-2i \sin(\pi k/n)} $$

$$ (44) \qquad = \frac{1}{4 \sin^2(\pi k/n)}. $$

Here we used the fact that the stationary distribution for the random walk on a ring is uniform [17], that is $\pi_q = 1/n$ for all $q$. Finally,

$$ (45) \qquad \gamma_k = \frac{1}{2 \sin(\pi k/n)} \underset{k \ll n}{\sim} \frac{n}{2 \pi k}. $$

In other word, the coupling coefficient scaling exponent is $\beta = -1$.

We observe that the continuous time random walk on a ring has the same scaling coefficients $(\alpha, \beta) = (2, -1)$ that the M/M/1 queue in the heavy traffic limit, and the PSD also scales as $1/f^{3/2}$. Intuitively this is because the eigenmodes of the semi-open string of the M/M/1 queue in the heavy traffic limit are circulating far away from the origin, just like the circulating counter-propagating eigenmodes on a ring (similarly to the circulating modes in a circular optical resonator). We confirmed this result by a numerical simulation Fig. 7.

<figure><img src="image_9.png" /><figcaption>FIGURE 7. Simulation of the continuous time random walk on a ring of size $n = 1000$. (a) Example of simulation trajectory (state on the ring versus time). (b) Periodogram (averaged over 30 realizations), the red line is the $1/f^{3/2}$ scaling.</figcaption></figure>

5.3. **Sequence of pulses:** $1/f^{1/2}$ noise. Existence of $1/f$ noise in pulse sequences has been demonstrated in [24], for example. Here we show that such signal can be obtained within the framework develop in this article. Until now we have only consider the case where the weight equals the state index ($x_i = i$), but other options are possible. By choosing a Kronecker delta for the graph weights: $x_q = \delta_q$, we obtain a pulsed signal where a pulse is emitted when the stochastic trajectory crosses the $n = 0$ state. Note that, because the number of internal states $n$ is larger than 2, this process is not a telegraph noise and, as we will see, the power spectral density is not Lorentzian.
The coupling coefficient simply satisfies:

$$ (46) \qquad \gamma_k^2 = \frac{1}{n^2} \sum_{q,r} \delta_q \delta_r w_n^{(r-q)k} = \frac{w_n^0}{n^2} = \frac{1}{n^2}. $$

Hence, there is no scaling of $\gamma_k$ with $k$, that is $\beta = 0$. This process has scaling exponents $(\alpha, \beta) = (2, 0)$, which satisfies $\alpha = 2 > |2\beta + 1| = 1$ and, by Criterion 6, the power spectral density scales as $S(\omega) \sim \omega^{-1/2}$, which is what we observe in the numerical simulation Fig. 8.

<figure><img src="image_10.png" /><figcaption>FIGURE 8. Simulation of a sequence of pulses with underlying random walk on a ring of size $n = 1000$. (a) Example of simulation trajectory. (b) Periodogram (averaged over 90 realizations), the green line is the $1/f^{1/2}$ scaling.</figcaption></figure>

## 6. CONCLUSION

Using a simple criterion relating the $1/f$ noise exponent to scalings of the eigenstructure of the infinitesimal generator, we observed that the M/M/1 queue in the heavy traffic and the random walk on a ring exhibit a $1/f^{3/2}$ noise. The existence of such a non-trivial behavior within the simplest model in queuing theory is an interesting example of the emergence of complex phenomena from simple rules. We also give an example of a sequence of pulses resulting from an underlying random walk on a ring and resulting in a $1/f^{1/2}$ noise.

The M/M/1 queue in the heavy traffic is a reflected Brownian motion [42], and interestingly $1/f$ noise has also been observed in other stochastic processes with limited or modified boundaries. We observed it for the random walk on a ring, that is a process with circular boundaries, but it has also been demonstrated in non-linear stochastic equations with reflective boundaries [33]. The reflective boundaries define the eigenmodes of the Markov chain graph (similarly to the eigenmodes of a resonator) and therefore the eigenstructure of the generator. The relationship between reflections in a Markov chain and $1/f$ noise is an interesting problem to explore.

As shown in App. F, the weight vector $\mathbf{x} = (x_i)$ can be understood as a signal on the transition graph, and the power spectral density is related to the graph Fourier transform of this signal. This shows that tools from the field of graph signal processing [31] can be used to study $1/f$ noise resulting from graph random walks. Conversely, the analysis of power spectral densities could be used to study graph signals. For example simulating a random walk on a large graph might be less computationally intensive than performing the full eigendecomposition and still provide useful insights by looking at the power spectral density and/or the autocorrelation function (power law scalings could help in classifying graphs).

## APPENDIX A. PROOF OF PROP. 1

*Proof.* Because the process is stationary $\langle X_{t+\tau} \rangle = \langle X_t \rangle$, and the autocorrelation function satisfies:

$$ (47) \qquad C_X(\tau) = \lim_{t \to \infty} [\langle X_{t+\tau} X_t \rangle - \langle X_t \rangle^2]. $$
Using the Markovian property:

$$ (48) \qquad P(X_{t+\tau} = j \wedge X_t = i) = P(X_{t+\tau} = j | X_t = i)P(X_t = i), $$

and because the CTMC is homogeneous

$$ (49) \qquad P(X_{t+\tau} = j | X_t = i) = P(X_{\tau} = j | X_0 = i) = \mathbf{P}_{ij}(\tau), $$

the first term is

$$ (50) \qquad \langle X_{t+\tau} X_t \rangle = \sum_{i,j} x_i x_j P(X_{t+\tau} = j \wedge X_t = i) = \sum_{i,j} x_i x_j \mathbf{P}_{ij}(\tau) P(X_t = i). $$

Moreover the chain is irreducible with a finite state space so it as an unique stationary distribution $\pi = (\pi_i)_i$ which satisfies $\lim_{t \to \infty} P(X_t = i) = \pi_i$, so that

$$ (51) \qquad \lim_{t \to \infty} \langle X_{t+\tau} X_t \rangle = \sum_{i,j} \pi_i x_i x_j \mathbf{P}_{ij}(\tau) = \langle \mathbf{x}, \mathbf{P}(\tau) \mathbf{x} \rangle_{\pi}. $$

The stationary average value of the process is

$$ (52) \qquad \lim_{t \to \infty} \langle X_t \rangle = \lim_{t \to \infty} \sum_i x_i P(X_t = i) = \sum_i x_i \pi_i = \langle \mathbf{x}, \mathbf{1} \rangle_{\pi}. $$

Since $\langle \mathbf{x}, \mathbf{1} \rangle^2 = \sum_{i,j} x_i x_j \pi_i \pi_j$, it results

$$ (53) \qquad C_X(\tau) = \sum_{i,j} x_i x_j \pi_i [\mathbf{P}_{ij}(\tau) - \pi_j]. $$

Moreover, $X$ is irreducible and has a stationary distribution, so the theorem of convergence to invariant distribution for a continuous time Markov chain is satisfied, and $\mathbf{P}_{ij}(\infty) \triangleq \lim_{\tau \to \infty} \mathbf{P}_{ij}(\tau) = \pi_j$. Therefore,

$$ (54) \qquad C_X(\tau) = \sum_{i,j} x_i x_j \pi_i [\mathbf{P}_{ij}(\tau) - \mathbf{P}_{ij}(\infty)] = \langle \mathbf{x}, [\mathbf{P}(\tau) - \mathbf{P}(\infty)] \mathbf{x} \rangle_{\pi}, $$

which is the expected result. $\square$

## APPENDIX B. PROOF OF PROP. 4

Here we explicit the relationship between the fundamental matrix and the infinitesimal generator. In this appendix $\sigma(\mathbf{A})$ designates the spectrum of the matrix $\mathbf{A}$, and $\mathbf{A} \oplus \mathbf{B}$ is the direct sum of $\mathbf{A}$ and $\mathbf{B}$. Moreover, $\mathbf{1}_1$ is the $1 \times 1$ identity matrix and $\mathbf{0}_1$ is the $1 \times 1$ zero matrix.

The result is straightforward to show if the generator is diagonalizable, as this is the case for a reversible process, however the result is valid for any generator of an irreducible chain. We start with the following lemma for the decomposition of the transition matrix, mostly resulting from the Perron-Frobenius theorem for transition matrices:

**Lemma 7** ([4], Lemma 8.2.2). If $\mathbf{P}$ is the transition matrix of an ergodic (irreducible) Markov chain then it exists an invertible matrix $\mathbf{S}$ so that $\mathbf{P} = \mathbf{S}[\mathbf{1}_1 \oplus \mathbf{K}]\mathbf{S}^{-1}$, with $1 \notin \sigma(\mathbf{K})$. Moreover, if the chain is regular then $\lim_{n \to \infty} \mathbf{K}^n = \mathbf{0}$.

As a corollary, we have an equivalent result for the core-nilpotent decomposition of the infinitesimal generator:

**Lemma 8.** If $\mathbf{G}$ is the generator of an irreducible HCTMC over a finite state space, then it exists an invertible matrix $\mathbf{S}$ so that $\mathbf{G} = \mathbf{S}[\mathbf{0}_1 \oplus \mathbf{\Gamma}]\mathbf{S}^{-1}$ where:

(i) $\mathbf{\Gamma}$ is invertible.

(ii) $\lim_{\tau \to \infty} e^{-\mathbf{\Gamma}\tau} = \mathbf{0}$.

(iii) for all $\tilde{\omega} \in \sigma(\mathbf{\Gamma})$, $\mathrm{Re}(\tilde{\omega}) > 0$.

*Proof.* Using Eq. (2) with $\tau = 1$, one has $\mathbf{P}(1) = e^{-\mathbf{G}}$. Since $\mathbf{P}(1)$ is a transition matrix of an irreducible chain, by Lemma 7 it exists an invertible matrix $\mathbf{S}$ so that $\mathbf{P}(1) = \mathbf{S}[\mathbf{1}_1 \oplus \mathbf{K}]\mathbf{S}^{-1}$, with $1 \notin \sigma(\mathbf{K})$.

Noting that $\mathbf{P}(1)$ is invertible (the inverse matrix being $e^{+\mathbf{G}}$), it follows that $\mathbf{K}$ is also invertible because $\det \mathbf{K} = \det \mathbf{P}(1) \neq 0$. Since the logarithm of an invertible matrix exists [16, Th. 1.27], but is not uniquely defined, let $-\mathbf{\Gamma}$ be a logarithm of $\mathbf{K}$ so that $\mathbf{K} = \exp(-\mathbf{\Gamma})$, and

$$ (55) \qquad \mathbf{P}(1) = \mathbf{S}[e^{-\mathbf{0}_1} \oplus e^{-\mathbf{\Gamma}}] \mathbf{S}^{-1} = \exp(-\mathbf{S}[\mathbf{0}_1 \oplus \mathbf{\Gamma}] \mathbf{S}^{-1}), $$
and therefore $\mathbf{G} = \mathbf{S}[\mathbf{0}_1 \oplus \mathbf{\Gamma}]\mathbf{S}^{-1}$. Moreover, since $1 \notin \sigma(\mathbf{K})$ then $0 \notin \sigma(\mathbf{\Gamma})$ and $\mathbf{\Gamma}$ is invertible, proving (i).

Given that a continuous time Markov chain over a finite state space is regular [19], Lemma 7 yields $\lim_{n \to \infty} \mathbf{K}^n = 0$ and (ii) follows

$$ (56) \qquad \lim_{\tau \to \infty} e^{-\Gamma\tau} = \lim_{n \to \infty} e^{-\Gamma n} = \lim_{n \to \infty} (e^{-\Gamma})^n = \lim_{n \to \infty} \mathbf{K}^n = \mathbf{0}. $$

Finally, to show (iii) we observe that if $\bar{\omega} \in \sigma(\mathbf{\Gamma})$ then $e^{-\bar{\omega}} \in \sigma(\mathbf{K})$ and $|e^{-\bar{\omega}}| < 1$, because $e^{-\bar{\omega}}$ is not the Perron-Frobenius eigenvalue of the stochastic matrix $\mathbf{P}(1)$, therefore $\mathrm{Re}(\bar{\omega}) > 0$. $\square$

We now use this result to prove Prop. 4:

*Proof.* In this proof, we make use of the group inverse $\mathbf{A}^{\sharp}$ of $\mathbf{A}$. It is a special case of the Drazin inverse for the matrix index of 1, which applies to us as we are dealing with matrices whose zero eigenvalue as multiplicity 1 (see [4] Theorem 8.2.1 and its proof). In other words, if $\mathbf{A}$ admits the decomposition $\mathbf{A} = \mathbf{S}[\mathbf{0}_1 \oplus \mathbf{Q}]\mathbf{S}^{-1}$ with $\mathbf{Q}$ invertible, then $\mathbf{A}^{\sharp} = \mathbf{S}[\mathbf{0}_1 \oplus \mathbf{Q}^{-1}]\mathbf{S}^{-1}$.

By Lemma 8, the transition matrix admits the following decomposition:

$$ (57) \qquad \mathbf{P}(\tau) = e^{-\mathbf{S}[\mathbf{0}_1 \oplus \mathbf{\Gamma}]\tau} = \mathbf{S}[\mathbf{1}_1 \oplus e^{-\mathbf{\Gamma}\tau}]\mathbf{S}^{-1}, $$

moreover using Lemma 8 (ii),

$$ (58) \qquad \mathbf{P}(\infty) = \lim_{\tau \to \infty} \mathbf{S}[\mathbf{1}_1 \oplus e^{-\mathbf{\Gamma}\tau}]\mathbf{S}^{-1} = \mathbf{S}[\mathbf{1}_1 \oplus \mathbf{0}]\mathbf{S}^{-1}, $$

so that $\mathbf{P}(\tau) - \mathbf{P}(\infty) = \mathbf{S}[\mathbf{0}_1 \oplus e^{-\mathbf{\Gamma}\tau}]\mathbf{S}^{-1}$. Therefore, the fundamental matrix is decomposed as

$$ (59) \qquad \mathbf{Z}(\omega) = \mathbf{S}[\mathbf{0}_1 \oplus \mathbf{J}(\omega)]\mathbf{S}^{-1}, $$

where

$$ (60) \qquad \mathbf{J}(\omega) \triangleq \int_0^\infty e^{-\mathbf{\Gamma}\tau} \cos(\omega\tau)d\tau = \frac{1}{2} \left[ \int_0^\infty e^{-\mathbf{\Gamma}\tau} e^{+i\omega\tau}d\tau + \int_0^\infty e^{-\mathbf{\Gamma}\tau} e^{-i\omega\tau}d\tau \right]. $$

Recall that the Laplace transform of the exponential of a matrix $\mathbf{X}$ is $\int_0^\infty e^{-ts}e^{t\mathbf{X}}dt = (s\mathbf{1} - \mathbf{X})^{-1}$, which is defined for all $s \in \mathbb{C} \setminus \sigma(\mathbf{X})$, and in particular

$$ (61) \qquad \int_0^\infty e^{-\mathbf{\Gamma}\tau} e^{\pm i\omega\tau} d\tau = (\mathbf{\Gamma} \mp i\omega\mathbf{1})^{-1}, $$

which is defined for all $\omega \in \mathbb{R}$ since by Lemma 8 (iii) $\sigma(\mathbf{\Gamma}) \not\subset i\mathbb{R}$. As a consequence, because $\mathbf{\Gamma}$ is invertible (Lemma 8 (i)),

$$ (62) \qquad \mathbf{J}(\omega) = \frac{1}{2} [(\mathbf{\Gamma} - i\omega\mathbf{1})^{-1} + (\mathbf{\Gamma} + i\omega\mathbf{1})^{-1}] = \mathbf{\Gamma} (\mathbf{\Gamma}^2 + \omega^2\mathbf{1})^{-1} = (\mathbf{\Gamma} + \omega^2\mathbf{\Gamma}^{-1})^{-1}. $$

Hence,

$$ (63) \qquad \mathbf{Z}(\omega) = \mathbf{S}[\mathbf{0}_1 \oplus (\mathbf{\Gamma} + \omega^2\mathbf{\Gamma}^{-1})^{-1}]\mathbf{S}^{-1} \triangleq \mathbf{M}^{\sharp}, $$

where $\mathbf{M}^{\sharp}$ is the group inverse of

$$ (64) \qquad \mathbf{M} = \mathbf{S}[\mathbf{0}_1 \oplus (\mathbf{\Gamma} + \omega^2\mathbf{\Gamma}^{-1})] \mathbf{S}^{-1} = \mathbf{S}[\mathbf{0}_1 \oplus \mathbf{\Gamma}] \mathbf{S}^{-1} + \omega^2 \mathbf{S}[\mathbf{0}_1 \oplus \mathbf{\Gamma}^{-1}] \mathbf{S}^{-1} = \mathbf{G} + \omega^2 \mathbf{G}^{\sharp}, $$

$$ \text{and finally } \mathbf{Z}(\omega) = (\mathbf{G} + \omega^2 \mathbf{G}^{\sharp})^{\sharp}. \qquad \square $$

## APPENDIX C. PROOF OF THEOREM. 5

*Proof.* Since the process is reversible, $\mathbf{G}$ is diagonalizable and has the spectral decomposition $\mathbf{G} = \sum_{k=1}^{n-1} \omega_k \mathbf{\Pi}_k$. Therefore, by Prop. 4, given that for all $\omega \in \mathbb{R}^+ \setminus \{0\}$ the function $x \mapsto x/(x^2 + \omega^2)$ is analytic on the domain $\mathbb{R}^+ \supset \sigma(\mathbf{G})$, the spectral decomposition of $\mathbf{Z}(\omega)$ is

$$ (65) \qquad \mathbf{Z}(\omega) = \sum_{k=1}^{n-1} \frac{\omega_k}{\omega_k^2 + \omega^2} \mathbf{\Pi}_k. $$

Because the process is reversible, it is stationary and applying Prop. 3 to the spectral decomposition of $\mathbf{Z}(\omega)$ provides Eq. (17).

The quadratic form $\langle \mathbf{x}, \mathbf{\Pi}_k \mathbf{x} \rangle_{\pi}$ is positive semi-definite because $\mathbf{\Pi}_k$ is a projector and $\sigma(\mathbf{\Pi}_k) = \{0, 1\} \subset \mathbb{R}^+$, so we can define $\gamma_k \triangleq \langle \mathbf{x}, \mathbf{\Pi}_k \mathbf{x} \rangle_{\pi}^{1/2}$. $\square$
APPENDIX D. PROOF OF CRITERION 6

We use Def. 3 of [37] as the definition for a function $f(x)$ to scale as $x^\alpha$. In the context we are interested in, it states that a function $f$ is asymptotically homogeneous at infinity if there exists $\alpha \in \mathbb{R}$ such that

$$ (66) \qquad \lim_{x \to \infty} \log \left( \frac{f(x)}{x^\alpha} \right) = 0. $$

The result derives from the asymptotic expansion for the power spectral density $S_X(\omega)$ (which we prove afterwards Sec. D.1):

**Proposition 9.** Assume the scalings $\omega_k \sim k^\alpha$ with $\alpha > 0$, and $\gamma_k \sim k^\beta$, so that $\alpha > |2\beta + 1|$, then

$$ (67) \qquad S_X(\omega) \sim \omega^\zeta - K\omega^{-2}, $$

where $\zeta = (2\beta - \alpha + 1)/\alpha$, and $K > 0$ if it exists $k \in \mathbb{N}$ so that $k = 2\beta + \alpha$, else $K = 0$.

Since $-2 < \zeta < 0$ then $\zeta + 2 > 0$ and it immediately follows:

$$ (68) \qquad \log \left( \frac{S_X(\omega)}{\omega^\zeta} \right) \sim \log \left( 1 - \frac{K}{\omega^{\zeta+2}} \right) \underset{\omega \to \infty}{\longrightarrow} 0. $$

Therefore $S_X$ is asymptotically homogeneous at $\omega \to \infty$ and $S_X(\omega)$ scales as $\omega^\zeta$.

D.1. **Proof of Prop. 9.**

*Proof.* We write the scaling laws as $\omega_k = \omega_0 k^\alpha$ and $\gamma_k = \gamma_0 k^\beta$, then $S_X(\omega) = \frac{\gamma_0^2}{\omega_0} S_{\alpha\beta}(\bar{\omega})$, where $\bar{\omega} \triangleq \omega/\omega_0$ and $S_{\alpha\beta}(\bar{\omega}) = \sum_{k=1}^\infty f_{\alpha\beta}(k)$, with

$$ (69) \qquad f_{\alpha\beta}(x) = \frac{x^{2\beta+\alpha}}{x^{2\alpha} + \bar{\omega}^2}. $$

We start with the Euler-Maclaurin asymptotic expansion to evaluate the sum:

$$ (70) \qquad \sum_{k=0}^{\infty} f_{\alpha\beta}(k) \sim I_{\alpha\beta}(\bar{\omega}) + \frac{f_{\alpha\beta}(\infty) + f_{\alpha\beta}(0)}{2} + \sum_{k=1}^{\infty} \frac{B_{2k}}{(2k)!} [f_{\alpha\beta}^{(2k-1)}(\infty) - f_{\alpha\beta}^{(2k-1)}(0)], $$

where $I_{\alpha\beta}(\bar{\omega}) \triangleq \int_0^\infty f_{\alpha\beta}(x) dx$, and $B_k$ is the $k$th Bernoulli number.

We first look at the behavior of $f_{\alpha\beta}$ and its derivatives as $x \to \infty$. Since $\alpha > 0$,

$$ (71) \qquad f_{\alpha\beta}(x) \underset{x \to \infty}{\sim} x^{2\beta-\alpha}, $$

and the $k$th derivative satisfies (see [20]):

$$ (72) \qquad f_{\alpha\beta}^{(k)}(x) \underset{x \to \infty}{\sim} \frac{\Gamma(2\beta + \alpha)}{\Gamma(2\beta + \alpha - k)} x^{2\beta - \alpha - k}, $$

where $\Gamma$ is the Euler gamma function. Because for all $k \in \mathbb{N}$, $2\beta - k < |2\beta + 1| < \alpha$, we have $f_{\alpha\beta}^{(k)}(\infty) \triangleq \lim_{x \to \infty} f_{\alpha\beta}^{(k)}(x) = 0$.

Moreover, because $\sum_{k=1}^\infty f_{\alpha\beta}(k) = \sum_{k=0}^\infty f_{\alpha\beta}(k) - f_{\alpha\beta}(0)$, the Euler-Maclaurin expansion Eq. (70) implies

$$ (73) \qquad S_{\alpha\beta}(\bar{\omega}) \sim I_{\alpha\beta}(\bar{\omega}) - R_{\alpha\beta}(\bar{\omega}), $$

where

$$ (74) \qquad R_{\alpha\beta}(\bar{\omega}) = \frac{1}{2} f_{\alpha\beta}(0) + \sum_{k=1}^{\infty} \frac{B_{2k}}{(2k)!} f_{\alpha\beta}^{(2k-1)}(0). $$

Therefore, by Lemma 10 and Lemma 11:

$$ (75) \qquad S_{\alpha\beta}(\bar{\omega}) = \frac{\pi}{2\alpha} \sec \left( \pi \frac{1+2\beta}{2\alpha} \right) \bar{\omega}^{\frac{2\beta-\alpha+1}{\alpha}} - K\bar{\omega}^{-2}. $$

□
Remark 5. The power spectral density $S_{\alpha\beta}(\bar{\omega})$ is finite at zero frequency ($\bar{\omega} = 0$), because

$$ (76) \qquad S_{\alpha\beta}(0) = \sum_{k=1}^{\infty} k^{2\beta-\alpha} = \zeta(\alpha - 2\beta), $$

where $\zeta$ is the Riemann zeta function. The series converges if $\alpha - 2\beta > 1$, which is satisfied because $\alpha > |2\beta + 1| \ge 2\beta + 1$.

**Lemma 10.** If $\alpha > 0$ and $\alpha > |2\beta + 1|$ then

$$ (77) \qquad I_{\alpha\beta}(\bar{\omega}) = \frac{\pi}{2\alpha} \sec\left(\pi \frac{1+2\beta}{2\alpha}\right) \bar{\omega}^{\frac{2\beta-\alpha+1}{\alpha}}, $$

where $\sec x \triangleq 1/\cos x$.

*Proof.* We use the fact that the fractional absolute moments of the Cauchy distribution:

$$ (78) \qquad p(x) = \frac{\gamma}{\pi} \frac{1}{x^2 + \gamma^2}, $$

are, for $|p| < 1$, $\langle|x|^p\rangle = \gamma^p \sec(\pi p/2)$. In other words,

$$ (79) \qquad \int_0^\infty \frac{x^p}{x^2 + \gamma^2} dx = \frac{\pi}{2} \gamma^{p-1} \sec\left(\frac{\pi p}{2}\right). $$

Because $\alpha > 0$, we can write:

$$ (80) \qquad I_{\alpha\beta}(\bar{\omega}) = \int_0^\infty \frac{x^{2\beta+\alpha}}{x^{2\alpha} + \bar{\omega}^2} dx = \int_0^\infty \frac{(x^\alpha)^{\frac{2\beta+\alpha+1}{\alpha}}}{(x^\alpha)^2 + \bar{\omega}^2} \frac{dx}{x}, $$

from the change of variable $y = x^\alpha$, we have $dy/y = \alpha dx/x$ and using Eq. (79) with $p = (2\beta+1)/\alpha$ it results

$$ (81) \qquad I_{\alpha\beta}(\bar{\omega}) = \frac{1}{\alpha} \int_0^\infty \frac{y^{\frac{2\beta+1}{\alpha}}}{y^2 + \bar{\omega}^2} dy = \frac{\pi}{2\alpha} \bar{\omega}^{\frac{2\beta+1}{\alpha} - 1} \sec\left(\pi \frac{1+2\beta}{2\alpha}\right), $$

which is defined because the condition $\alpha > |2\beta + 1|$ implies $|p| < 1$. $\square$

Remark 6. This proof exhibits a relationship between the power spectral density of a Markov chain with power law scaling eigenstructure and the fractional absolute moments of the Cauchy distribution, basically $S(\omega) \sim \langle|x|^{1+\zeta}\rangle \sim \omega^\zeta$.

Note that the fractional absolute moments are closely related to the Mellin transform which is often "Used in place of Fourier's transform when scale invariance is more relevant than shift invariance" [2]. It would therefore be interesting to study the potential relation between $1/f$ power spectral densities and the Mellin transform. Maybe one could derive a variant of the Wiener-Kinchine theorem based on the Mellin transform instead of the Fourier transform?

**Lemma 11.** If $\alpha > 0$ then $R_{\alpha\beta}(\bar{\omega}) = K\bar{\omega}^{-2}$, where

$$ (82) \qquad K = \begin{cases} \frac{1}{2} & \text{if } 2\beta + \alpha = 0, \\ \frac{B_{2p}}{(2p)!} \Gamma(2\beta + \alpha) & \text{if it exists } p \in \mathbb{N} \setminus \{0\} \text{ so that } 2p = 2\beta + \alpha + 1, \\ 0 & \text{otherwise.} \end{cases} $$

That is $K > 0$ if it exists $k \in \mathbb{N}$ so that $k = 2\beta + \alpha$, else $K = 0$.

*Proof.* Since $\alpha > 0$,

$$ (83) \qquad f_{\alpha\beta}(x) \underset{x\to 0}{\sim} \frac{x^{2\beta+\alpha}}{\bar{\omega}^2}, $$

and $k$th derivative ($k \in \mathbb{N}$) satisfies:

$$ (84) \qquad f_{\alpha\beta}^{(k)}(x) \underset{x\to 0}{\sim} \frac{\Gamma(2\beta + \alpha) x^{2\beta + \alpha - k}}{\Gamma(2\beta + \alpha - k) \bar{\omega}^2}. $$

Therefore, for all $k \in \mathbb{N}$,

$$ (85) \qquad f_{\alpha\beta}^{(k)}(0) = \begin{cases} \Gamma(2\beta + \alpha)/\bar{\omega}^2 & \text{if } k = 2\beta + \alpha, \\ 0 & \text{otherwise.} \end{cases} $$
As a consequence,

$$ (86) \qquad R_{\alpha\beta}(\bar{\omega}) = \begin{cases} \frac{1}{2\bar{\omega}^2} & \text{if } 2\beta + \alpha = 0, \\ \frac{B_{2p} \Gamma(2\beta+\alpha)}{(2p)! \bar{\omega}^2} & \text{if it exists } p \in \mathbb{N} \setminus \{0\} \text{ so that } 2p = 2\beta + \alpha + 1, \\ 0 & \text{otherwise.} \end{cases} $$

□

## APPENDIX E. PSD OF BIRTH-DEATH PROCESSES IN THE HEAVY TRAFFIC REGIME

A birth-death process is a Markov process where from a given state one can only realize a birth transition increasing the state count, or a death transition reducing the state count (except from the first state from which on ly births are allowed). The infinitesimal generator of this process is a tridiagonal matrix:

$$ (87) \qquad \mathbf{G} = \begin{pmatrix} \lambda_0 & -\lambda_0 & & & \\ -\mu_1 & \mu_1 + \lambda_1 & -\lambda_1 & & \\ & -\mu_2 & \mu_2 + \lambda_2 & -\lambda_2 & \\ & & \ddots & \ddots & \\ & & & \ddots & \ddots \end{pmatrix}, $$

where $\lambda_n$ are the birth rates and $\mu_n$ the death rates. In this appendix, we provide a general formula for the power spectral density of a birth-death process in the heavy traffic regime. More precisely, we see that the coupling coefficient is related to the family of orthogonal polynomials associated to the characteristic polynomial of the generator.

**E.1. Characteristic polynomial of the infinitesimal generator.** Because the matrix is tridiagonal, the characteristic polynomial $f_n$ of the truncated $n \times n$ generator $G_n$ satisfies the three terms recurrence relation:

$$ (88) \qquad x f_n(x) = -\lambda_n \mu_n f_{n-1}(x) + (\mu_{n+1} + \lambda_{n+1}) f_n(x) - f_{n+1}(x), $$

with initial conditions $f_0(x) = 1$ and $f_1(x) = \lambda_1 - x$. Since $\lambda_n \mu_n > 0$, from Favard theorem, the polynomial $f_n$ is orthogonal [25]. Therefore the eigenvalues of $G$ are roots of orthogonal polynomials.

Moreover the roots of successive polynomials $f_n$ and $f_{n+1}$ from a family of orthogonal polynomials $\{f_n\}_{n \in \mathbb{N}}$, which are the eigenvalues of $G_n$ and $G_{n+1}$, are interleaved [25]. In the limit $n \to \infty$, the eigenfrequencies of $G$ are therefore densely covering a wide span of values, which is favorable to the emergence of $1/f$ scaling.

**E.2. Heavy traffic limit.** We define the server utilization $\rho_i \triangleq \lambda_{i-1}/\mu_i$. In the heavy traffic approximation $\rho_i \sim 1$, for all $i$, and $G$ is symmetric.

In [29] it is shown that the coefficients of the eigenvector of a tridiagonal matrix are related to the family of orthogonal characteristic polynomials $\{f_n\}_{1 \le i \le n}$. More specifically, for symmetric tridiagonal matrix the eigenvector $v_k = (v_{k,i})_{1 \le i \le n}$ associated to the eigenvalue $\omega_k$ is:

$$ (89) \qquad v_{k,i} = q_k \psi_i(\omega_k), $$

where $\psi_i(x) \triangleq f_{i-1}(x)/\beta_{i-1}$ and $\beta_i \triangleq \prod_{l=1}^i \mu_l$. The normalization coefficient is:

$$ (90) \qquad q_k^2 = \frac{\beta_{n-1}}{\psi_n(\omega_k) \prod_{i \neq k} (\omega_k - \omega_i)}. $$

The coupling factor is therefore:

$$ (91) \qquad \gamma_k^2 = q_k^2 \left( \sum_{i=1}^{n} \pi_i x_i \psi_i(\omega_k) \right) \left( \sum_{i=1}^{n} x_i \psi_i(\omega_k) \right), $$

The steady-state population for the state $i$ of the birth-death process is $\pi_i = \pi_0 \prod_{l=1}^k \rho_i$, so in the heavy traffic limit the populations are almost equidistributed ($\pi_i \sim \pi_0$ for all $i$). Therefore,

$$ (92) \qquad \gamma_k \approx \sqrt{\pi_0 q_k \hat{x}(\omega_k)}, $$
where we introduced the generalized Fourier serie with respect to the family of functions {$\psi_i$}$_{1\le i\le n}$:

$$
(93) \qquad \hat{\mathbf{x}}(\omega) \triangleq \sum_{i=1}^{n} x_i \psi_i(\omega).
$$

We see in App. F that, when the generator is symmetric, the factor $\gamma_k$ is always a Fourier serie,
more specifically it is the graph Fourier transform of the Markov chain graph.

APPENDIX F. RELATIONSHIP BETWEEN THE PSD AND THE GRAPH FOURIER TRANSFORM

A graph can be associated to the Markov where each vertice is a state of the chain and the
edges of the graph are defined by the transition rates. When the generator **G** is symmetric, it is
also the Laplacian matrix of the graph associated to the Markov chain. Let $\omega_k$ be an eigenvalue
of **G** and $\mathbf{v}_k = (v_{k,i})_{1\le i\le n}$ the corresponding eigenvector, the *graph Fourier transform* [31] $\hat{\mathbf{x}}$ of
$\mathbf{x} \in \mathbb{R}^n$ at $\omega_k$ is defined as

$$
(94) \qquad \hat{\mathbf{x}}(\omega_k) \triangleq \langle \mathbf{x}, \mathbf{v}_k \rangle = \sum_{i=1}^{n} x_i v_{k,i}.
$$

The coupling factor $\gamma_k$ satisfies

$$
(95) \qquad \gamma_k^2 = \langle \mathbf{x}, \mathbf{v}_k^\top \mathbf{v}_k \mathbf{x} \rangle_\pi = \left( \sum_{i=1}^n \pi_i x_i v_{k,i} \right) \left( \sum_{i=1}^n x_i v_{k,i} \right).
$$

Since the generator is symmetric, the transition matrix **P** is a doubly stochastic matrix and the Markov process as a uniform limiting distribution: $\pi_i = 1/n$ for all $i$. Therefore, the coupling factor can be written using the graph Fourier transform Eq. (94):

$$
(96) \qquad \gamma_k = \frac{1}{\sqrt{n}} \sum_{i=1}^{n} x_i v_{k,i} = \frac{1}{\sqrt{n}} \hat{\mathbf{x}}(\omega_k).
$$

We can define the *graph power spectral density* as

$$
(97) \qquad S_X(\omega) = \frac{1}{n} \sum_{k=1}^{n} \hat{\mathbf{x}}(\omega_k)^2 \frac{\omega_k}{\omega_k^2 + \omega^2}.
$$

REFERENCES

1. P. Berndt and D. T. Kaplan, *Nonstationarity and 1/f noise characteristics in heart rate*, Am. J. Physiol. - Regulatory Integrative Comp. Physiol. **276** (1999), R1–R9.

2. J. Bertrand, P. Bertrand, and J.-P. Ovarlez, *The mellin transform*, 1995.

3. A. E. Brouwer and W. H. Haemers, *Spectra of graphs*, Springer New York, NY, 2012.

4. S. L. Campbell and C. D. Meyer, *Generalized inversed of linear transformations*, Pitman, London, 1979.

5. P. Coolen-Schrijner and E. A. van Doorn, *The deviation matrix of a continuous-time markov chain*, Probab. Eng. Inf. Sci. **16** (2002), 351–366.

6. I. Csabai, *1/f noise in computer network traffic*, J. Phys. A: Math. Gen. **27** (1994), L417–L421.

7. J. Dalibard and C. Cohen-Tannoudji, *Dressed-atom approach to atomic motion in laser light: the dipole force revisited*, J. Opt. Soc. Am. B **2** (1985), 1707.

8. S. Erland and P. E. Greenwood, *Constructing $1/\omega^\alpha$ noise from reversible markov chains*, Phys. Rev. E **76** (2007), 031114.

9. S. Erland, P. E. Greenwood, and L. M. Ward, '1/f<sup>α</sup> noise' is equivalent to an eigenstructure power relation, Europhysics Letters **95** (2011), 60006.

10. T. Geisel, A. Zacherl, and G. Radons, *Chaotic diffusion and 1/f-noise of particles in two-dimensional solids*, Z. Physik B - Condensed Matter **71** (1988), 117–127.

11. D. L. Gilden, T. Thornton, and M. W. Mallon, *1/f noise in human cognition*, Science **267** (1995), 1837–1839.

12. D. T. Gillespie, *A general method for numerically simulating the stochastic time evolution of coupled chemical reactions*, Journal of Computational Physics **22** (1976), 403–434.

13. A. Greco, S. Perri, S. Servidio, E. Yordanova, and P. Veltri, *The complex structure of magnetic field discontinuities in the turbulent solar wind*, The Astrophysical Journal Letters **823** (2016), L39.

14. L. Z. Hadid, F. Sahraoui, K. H. Kiyani, A. Retinò, R. Modolo, P. Canu, A. Masters, and M. K. Dougherty, *Nature of the mhd and kinetic scale turbulence in the magnetosheath of saturn: Cassini observations*, The Astrophysical Journal Letters **813** (2015), L29.

15. M. Haviv, *Queues*, Springer New York, NY, 2013.

16. N. J. Higham, *Functions of matrices*, Society for Industrial and Applied Mathematics, 2008.

17. Marcus M (https://math.stackexchange.com/users/215322/marcus m), *Stationary distribution and limiting distribution of a random walk*, Mathematics Stack Exchange, URL:https://math stackexchange com/q/2947045 (version: 2018-10-08).
18. D. Fulea (https://math stackexchange.com/users/550003/dan fulea), If p is the transition matrix of a reversible markov chain, why are its eigenvalues real?, Mathematics Stack Exchange, URL:https://math stackexchange.com/q/3157533 (version: 2019-03-21).

19. Arthur (https://math stackexchange.com/users/567412/arthur), Regular continuous-time markov chain with finite state space, Mathematics Stack Exchange, URL:https://math stackexchange.com/q/4025628 (version: 2021-02-14).

20. Clement C. (https://math stackexchange.com/users/75808/clement c), Differentiation and asymptotic equivalents, Mathematics Stack Exchange, URL:https://math stackexchange.com/q/2381258 (version: 2020-06-12).

21. D. L. Iglehart and W. Whitt, Multiple channel queues in heavy traffic. i, Adv. Appl. Prob. **2** (1970), 150–177.

22. D.-Q. Jiang and F.-X. Zhang, The green-kubo formula and power spectrum of reversible markov processes, J. Math. Phys. **44** (2003), 4681–4689.

23. J. F. C. Kingman, On queues in heavy traffic, Journal of the Royal Statistical Society: Series B (Methodological) **24** (1962), 383–392.

24. A. Kononovicius and B. Kaulakys, 1/f noise from the sequence of nonoverlapping rectangular pulses, Phys. Rev. E **107** (2023), 034117.

25. T. H. Koornwinder, Orthogonal polynomials, pp. 145–170, Springer Vienna, 2013.

26. C. Leyris, F. Martinez, M. Valenza, A. Hoffmann, J. C. Vildeuil, and F. Roy, Impact of random telegraph signal in cmos image sensors for low-light levels, oct 2006, pp. 376 – 379.

27. L. Mandel and E. Wolf, Optical coherence and quantum optics, Cambridge University Press, 1995.

28. E. Milotti, 1/f noise: a pedagogical review, apr 2002, arXiv:physics/0204033.

29. L. G. Molinari, An introduction to random matrices: Tridiagonal matrices, http://pcteserver.mi.infn.it/~molinari/RMT/RM_tridiagonal.pdf, 2019, Accessed: 2022-08-10.

30. S. Noschese, L. Pasquini, and L. Reichel, Tridiagonal toeplitz matrices: Properties and novel applications, Numer. Linear Algebra Appl. **20** (2013), 302–326.

31. B. Ricaud, P. Borgnat, N. Tremblay, P. Gonçalves, and P. Vandergheynst, Fourier could be a data scientist: From graph fourier transform to signal processing on graphs, Comptes Rendus Physique **20** (2019), 474–488.

32. J. Ruseckas and B. Kaulakys, 1/f noise from nonlinear stochastic differential equations, Phys. Rev. E **81** (2010), 031105.

33. —, Scaling properties of signals as origin of 1/f noise, J. Stat. Mech. **2014** (2014), P06005.

34. M. Takayasu, Dynamic complexity in the internet traffic, pp. 329–358, Springer Berlin Heidelberg, 2005.

35. A. van der Ziel, Unified presentation of 1/f noise in electron devices: fundamental 1/f noise sources, Proceedings of the IEEE **76** (1988), 233–258.

36. T. Vanderbruggen and M. W. Mitchell, Near-resonant optical forces beyond the two-level approximation for a continuous source of spin-polarized cold atoms, Phys. Rev. A **87** (2013), 033410.

37. F. Vigneron, On the wiener-khinchin transform of functions that behave as approximate power-laws. applications to fluid turbulence., 2019, arXiv:1909.06078.

38. R. F. Voss and J. Clarke, '1/f noise' in music: Music from 1/f noise, J. Acoust. Soc. Am **63** (1978), 258.

39. A. A. Walma, On the superposition of relaxation spectra as an explanation for 1/f noise, Revue Phys. Appl. **15** (1980), 1435–1443.

40. L. M. Ward and P. E. Greenwood, 1/f noise, Scholarpedia **2** (2007), 1537.

41. W. Whitt, Asymptotic formulas for markov processes with applications to simulation, Operations Research **40** (1992), 199–415.

42. Wikipedia contributors, Reflected brownian motion — Wikipedia, the free encyclopedia, 2022, [Online; accessed 5-January-2023].

43. A. C. Yadav and N. Kumar, Scaling theory for the 1/f noise, mar 2021, arXiv:2103.11608.

44. Y. Yamamoto, Fundamentals of noise processes: 1/f noise and random telegraph signals, https://www.nii.ac.jp/qis/first-quantum/e/forStudents/lecture/pdf/noise/chapter9.pdf, Accessed: 2022-08-23.

45. D. D. Yao, First-passage-time moments of markov processes, Journal of Applied Probability **22** (1985), 939–945.

KOHERON
3 RUE DU SOUS MARIN VENUS
56100 LORIENT
FRANCE

Email address: thomas.vanderbruggen@koheron.com