Onrstein-Uhlenbeck Process with power-law-correlated noise It's a copy I posted as: https://math.stackexchange.com/questions/4356212/orstein-uhlenbeck-process-with-power-law-correlated-noise
Consider a noise-driven drifting system given by the Langevin Eq:
$$\dot{v}+\gamma v =\xi$$
Where $v$ denotes velocity with $\gamma$ being drift coefficient and $\xi$ a stochastic-noisy drive, but unlike Orstein-Uhlenbeck Process, $\xi 's$ are not $uncorrelated$ rather they themselves are $\frac{1}{f}$ noise with power-law autocorrelation function [recalling from Wiener-Khinchin theorem that the autocorrelation function is fourier transform of the power spectral density where $\frac{1}{f}$ noise has Power Spectral density $S(f)\propto \frac{1}{f^{\alpha}}$ with $1 \leq \alpha \leq 2$], say (assuming stationarity):
$$\langle \xi(t) \ \xi(0) \rangle = \tilde{\Gamma}|t|^{-\theta} \ \ ; \ \ \ 0 \leq \theta \leq 1$$
Therefore the velocity profile of $v(t)$ is as:
$$v(t)= \Big[ v_0 +  \int^{t}_{0}e^{\gamma t'} \xi(t')dt' \Big] e^{-\gamma t} $$ where $\langle \xi(t) \rangle = 0$ will be considered, thus the average $\langle v(t) \rangle$ is same as Orstein-Uhlenbeck case,
BUT, what about the case of average energy, or second moment $\langle v^2(t) \rangle$, i.e.:
$$\langle v^2(t) \rangle = v_0^2 e^{-2 \gamma t} + e^{- 2 \gamma t } \int^{t}_{0} dt_1 \int^{t}_{0} dt_2 e^{\gamma{(t_1 + t_2)}} \langle \xi(t_1) \xi(t_2) \rangle $$
This expression is well treatable with delta correlated gaussian white noise, but it's nontrivial here.
How do I get passed with this integral ?
$$\int^{t}_{0} dt_1 \int^{t}_{0} dt_2 e^{\gamma{(t_1 + t_2)}}|t_1 - t_2|^{-\theta} $$
I am stuck here. Kindly tell me how to handle such long-range-correlated-noisy setup.
 A: Let $S_{\xi}(u)$ denote the spectrum of $\xi (t)$, that is if $K_{\xi}(t'-t'')=\langle \xi(t') \xi(t'')\rangle$ then $S_{\xi}(u)=\int_{-\infty}^{\infty}dt e^{-\mathfrak j 2\pi u t}K_{\xi}(t)$ and $K_{\xi}(t)=\int_{-\infty}^{\infty}du e^{\mathfrak j 2\pi u t}S_{\xi}(u)$.
In the integral  $$M(\theta) = \int^{t}_{0} dt_1 \int^{t}_{0} dt_2 e^{\gamma{(t_1 + t_2)}}|t_1 - t_2|^{-\theta} \tag{1}$$ replace the autocorrelation kernel $K_{\xi}(t)=|t|^{-\theta}$with its Fourier transform, ie. its spectrum and interchange  the order of integration. According to Item 311 in wiki the Fourier transform of $|x|^{\alpha}$, $-1<\alpha <0$ is
$-\frac{2 sin(\pi\alpha/2)\Gamma(1+\alpha)}{|2\pi \xi|^{1+\alpha} }$ which in your notation use $-\theta \leftarrow  \alpha$, $u \leftarrow  \xi$ snd $t \leftarrow  x$, therefore $$S_{\xi}=\frac{2 sin(\pi\theta/2)\Gamma(1-\theta)}{|2\pi u|^{1 -\theta} } \tag{S}\label{S}$$
$$M(\theta) = \int^{t}_{0} dt_1 \int^{t}_{0} dt_2 e^{\gamma{(t_1 + t_2)}}|t_1 - t_2|^{-\theta}$$ \
$$= \int^{t}_{0} dt_1 \int^{t}_{0} dt_2 e^{\gamma{(t_1 + t_2)}}\int_{-\infty}^{\infty}du e^{\mathfrak j 2\pi u(t_1-t_2)}S_{\xi}(u)\\
=\int_{-\infty}^{\infty}du S_{\xi}(u) \int^{t}_{0} dt_1 \int^{t}_{0} dt_2 e^{\gamma{(t_1 + t_2)}}e^{\mathfrak j 2\pi u (t_1-t_2)}\\
=\int_{-\infty}^{\infty} du S_{\xi}(u) \int^{t}_{0} dt_1 e^{\gamma{t_1}} e^{\mathfrak j 2\pi u t_1} \int^{t}_{0} dt_2 e^{\gamma t_2} e^{-\mathfrak j 2\pi u t_2} $$
Now $\int^{t}_{0} dt' e^{\gamma t'} e^{\mathfrak j 2\pi u t'}=-\frac{1-e^{(\gamma +\mathfrak j 2\pi u )t}}{\gamma+\mathfrak j 2\pi u}$, therefore
$$p(u,t)=\int^{t}_{0} dt_1 e^{\gamma{t_1}} e^{\mathfrak j 2\pi u t_1} \int^{t}_{0} dt_2 e^{\gamma t_2} e^{-\mathfrak j 2\pi u t_2}= \frac{1-e^{(\gamma +\mathfrak j 2\pi u) t}}{\gamma+\mathfrak j 2\pi u}\frac{1-e^{(\gamma -\mathfrak j 2\pi u )t}}{\gamma-\mathfrak j 2\pi u}\\
=\left|\frac {1-e^{(\gamma +\mathfrak j 2\pi u) t}}{\gamma+\mathfrak j 2\pi u}\right|^2\\
=\left|\frac{1+e^{2\gamma t} -2e^{\gamma t} cos(2\pi ut)}{\gamma^2+(2\pi u)^2}\right|^2 \tag{p}\label{p}$$ and we get with $\eqref{S}$ and $\eqref{p}$
$$M(\theta) = \int_{-\infty}^{\infty} du S_{\xi}(u)p(u,t)\\
=\int_{-\infty}^{\infty} du \frac{2 sin(\pi\theta/2)\Gamma(1-\theta)}{|2\pi u|^{1 -\theta} }\left|\frac{1+e^{2\gamma t} -2e^{\gamma t} cos(2\pi ut)}{\gamma^2+(2\pi u)^2}\right|^2\tag{M}\label{M}$$
Since the integral $\int_{-A}^A du \frac{1}{|u|^{1 -\theta}}=2\int_{0}^A du (u^{-1 +\theta})=\frac{A^\theta}{\theta}$ is finite for all finite $A$ upper limit, and therefore the integral in the variable $u$ in $\eqref{M}$ is finite, convergent, for all $t$ and can be calculated numerically.
