For a particle in a harmonic potential with a damping rate due to the enviroment $\Gamma$, a natural frequency $\omega_0$ and a corresponding trap stiffness $k_0 = m\omega_0^2$ with a mass $m$ and an external force acting on the particle $F_{th}$ the equation of motion for 1 degree of freedom is:
$$ m\ddot{x}(t) + m\Gamma \dot{x}(t) + k_0 x(t) = F_{th}(t)$$replacing $k_0$ with $m\omega^2$ we get
$$ \ddot{x}(t) + \Gamma \dot{x}(t) + \omega_0^2 x(t) = \dfrac{F_{th}(t)}{m}$$If we then use the Fourier relationship
$$ \dfrac{d^n f(t)}{dt^n} = (i\omega)^n \tilde{f}(\omega) $$Then the Fourier transform of the equation of motion is:
$$ (i\omega)^2 \tilde{x}(\omega) + \Gamma_0(i\omega)\tilde{x}(\omega) + \omega_0^2 \tilde{x}(\omega) = \dfrac{\tilde{F_{th}}(\omega)}{m} $$which is equal to:
$$ -\omega^2 \tilde{x}(\omega) + \Gamma_0(i\omega)\tilde{x}(\omega) + \omega_0^2 \tilde{x}(\omega) = \dfrac{\tilde{F_{th}}(\omega)}{m} $$Rearranging this for $\tilde{x}(\omega)$ we get:
$$ \tilde{x}(\omega) = \dfrac{1}{m} \dfrac{\tilde{F_{th}}(\omega)}{\omega_0^2-\omega^2 + i(\Gamma_0\omega)} $$The formula for the power spectral density (the Fourier transform of the autocorrelation function) is:
$$ S_{xx}(\omega) = \int_{-\infty}^{+\infty} \langle x(t) x(t+\tau)^* \rangle e^{-i\omega \tau} d\tau $$From https://journals.aps.org/rmp/pdf/10.1103/RevModPhys.86.1391 using the Wiener-Khinchin theorem.
And the formula for the reverse Fourier transform is
$$ x(t) = \dfrac{1}{\sqrt{2\pi}} \int_{-\infty}^{+\infty} \tilde{x}(\omega)e^{+i\omega t} d\omega $$Now using $\dfrac{1}{2\pi} \int_{-\infty}^{\infty} e^{-i\omega\tau} e^{i\omega_2\tau} d\tau = \dfrac{1}{2\pi} \int_{-\infty}^{\infty} e^{i(\omega_2-\omega)\tau} d\tau = \delta(\omega_2 - \omega)$ from the definition of the delta function
If we now use the result, which is derived below here, that:
$$ \langle \tilde{F_{th}}(\omega_1) \tilde{F_{th}}(\omega) \rangle = \sqrt{\dfrac{\pi}{2}} \delta(\omega_1 + \omega) \tilde{f}(\dfrac{\omega_1 - \omega}{2}) $$Where $\tilde{f}$ is the Fourier transform of the 2 point correlation function of the thermal noise $F_{th}$.
We get:
$$ S_{xx}(\omega) = \dfrac{1}{m^2} \int_{-\infty}^{\infty} d\omega_1 e^{+i(\omega + \omega_1) t} \dfrac{\delta(\omega_1 + \omega) \sqrt{\dfrac{\pi}{2}} \tilde{f}(\frac{\omega_1 - \omega}{2})}{(\omega_0^2-\omega_1^2 + i(\Gamma_0\omega_1))\cdot(\omega_0^2-\omega^2 - i(\Gamma_0\omega))} $$Using this result we get (using $\omega_1 = - \omega$ for each term $\omega_1$ inside the integral - due to the delta function):
$$ S_{xx}(\omega) = \dfrac{1}{m^2} \dfrac{\sqrt{\dfrac{\pi}{2}} \tilde{f}(-\omega)}{(\omega_0^2-\omega^2 - i(\Gamma_0\omega))\cdot(\omega_0^2-\omega^2 + i(\Gamma_0\omega))} $$which when expanded comes to:
$$ S_{xx}(\omega) = \sqrt{\dfrac{\pi}{2}} \dfrac{1}{m^2} \dfrac{\tilde{f}(-\omega)}{\omega_0^4 - 2\omega_0^2 \omega^2 + \omega^4 + \Gamma_0^2\omega^2} = \sqrt{\dfrac{\pi}{2}} \dfrac{1}{m^2} \dfrac{\tilde{f}(-\omega)}{(\omega_0^2 - \omega^2)^2 + (\Gamma_0\omega)^2} $$$\tilde{f}$, as mentioned earilier, is the Fourier transform of the 2 point correlation function of the thermal noise $F_{th}$. Taking this to be constant is equivalent to assuming that the noise is Markovian i.e. that it is uncorrelated from time $t$ to time $t + \tau$.
This is the case because the reverse Fourier transform of constant $\tilde{f}$ results in $f(s) = \dfrac{1}{\sqrt{2\pi}}\int{\tilde{f}e^{-i\omega s} d\omega} = \tilde{f}\delta(s)$ which is the two-point correlation function in time for Markovian noise. That is to say, it is uncorrelated in time.
therefore we get the result that $$ S_{xx}(\omega) = \sqrt{\dfrac{\pi}{2}} \dfrac{1}{m^2} \dfrac{\tilde{f}}{(\omega_0^2-\omega^2)^2 + (\Gamma_0\omega)^2} $$
Using the Fourier transform $\tilde{h}(\omega) = \dfrac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} x(t) e^{-i\omega t} dt$
Defining $$ \langle h(t_1)h(t_2) \rangle = f(t_1 - t_2) $$ which is the two-point correlation function of the thermal noise $F_{th}$.
we get
$$ = \dfrac{1}{2\pi} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} dt_1 dt_2 f(t_1-t_2) e^{-i\omega_1 t_1} e^{-i\omega_2 t_2} $$Now defining $t_+ = t_1 + t_2$ and $t_- = t_1 - t_2$ and replacing $t_1$ with $\dfrac{1}{2}(t_+ + t_-)$ and $t_2$ with $\dfrac{1}{2}(t_+ - t_-)$ in the integral we get:
$$ = \dfrac{1}{4\pi} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} dt_+ dt_- f(t_-) e^{-\frac{1}{2}i\omega_1 (t_+ + t_-)} e^{-\frac{1}{2}i\omega_2 (t_+ - t_-)} $$By taking the first term to be $$\dfrac{1}{2\pi} \int_{-\infty}^{\infty} dt_+ e^{i t_+ (\frac{-\omega_1 - \omega_2}{2})} = \delta(\dfrac{-\omega_1 - \omega_2}{2}) $$ from the definition of the delta function
And the second term to be $$\int_{-\infty}^{\infty} dt_- e^{-i t_- (\frac{\omega_1 - \omega_2}{2})} f(t_-) = \sqrt{2\pi} \tilde{f}(\dfrac{\omega_1-\omega_2}{2}) $$ from the definition of the Fourier transform
$\delta(\dfrac{-\omega_1 - \omega_2}{2})$ is equivalent to $\delta(\omega_1 + \omega_2)$ which leaves us with:
The form of the constant $\tilde{f}$ can be derived readily by realising that the integral of a power spectral density yields the positional variance of the oscillator, namely $\langle x^2 \rangle$. As used in https://link.aps.org/doi/10.1103/RevModPhys.86.1391
Therefore
$$ \langle x^2 \rangle = \int_0^{\infty} S_{xx}(\omega) \dfrac{d\omega}{2\pi} $$$$ = \sqrt{\dfrac{\pi}{2}} \dfrac{\tilde{f}}{m^2} \dfrac{1}{2\pi} \int_{0}^{\infty} \dfrac{1}{(\omega_0^2-\omega^2)^2 + (\Gamma_0\omega)^2} d\omega $$$$ = \sqrt{\dfrac{1}{\sqrt{2\pi}}} \dfrac{\tilde{f}}{2m^2} \int_{0}^{\infty} \dfrac{1}{(\omega_0^2-\omega^2)^2 + (\Gamma_0\omega)^2} d\omega $$This integral can be calculated like so:
Preliminaries
Using $a^2+b^2=(a+ib)(a-ib)$, we get
\begin{align} \left(\omega_0^2-\omega^2\right)^2+(\Gamma_0\omega)^2 &=\left(\omega^2+i\Gamma_0\omega-\omega_0^2\right)\left(\omega^2-i\Gamma_0\omega-\omega_0^2\right) \end{align}The roots for the two quadratics are $$ \omega=\frac{\pm i\Gamma_0\pm\sqrt{4\omega_0^2-\Gamma_0^2}}2 $$ Thus, re-pairing the roots gives
\begin{align} &\left(\omega_0^2-\omega^2\right)^2+(\Gamma_0\omega)^2\\ &=\left(\omega^2+\omega\sqrt{4\omega_0^2-\Gamma_0^2}+\omega_0^2\right)\left(\omega^2-\omega\sqrt{4\omega_0^2-\Gamma_0^2}+\omega_0^2\right)\\ &=\left[\left(\omega+\sqrt{\omega_0^2-\tfrac14\Gamma_0^2}\right)^2+\tfrac14\Gamma_0^2\right]\left[\left(\omega-\sqrt{\omega_0^2-\tfrac14\Gamma_0^2}\right)^2+\tfrac14\Gamma_0^2\right] \end{align}Then Partial Fractions gives
\begin{align} &\frac1{\left(\omega_0^2-\omega^2\right)^2+(\Gamma_0\omega)^2}\\[12pt] &=\frac1{2\omega_0^2}\left[\frac{1+\frac\omega{\sqrt{4\omega_0^2-\Gamma_0^2}}}{\left(\omega+\sqrt{\omega_0^2-\tfrac14\Gamma_0^2}\right)^2+\tfrac14\Gamma_0^2} +\frac{1-\frac\omega{\sqrt{4\omega_0^2-\Gamma_0^2}}}{\left(\omega-\sqrt{\omega_0^2-\tfrac14\Gamma_0^2}\right)^2+\tfrac14\Gamma_0^2}\right] \end{align}Integration
Due to the evenness of the integrand,
\begin{align} &\int_0^\infty\frac1{\left(\omega_0^2-\omega^2\right)^2+(\Gamma_0\omega)^2}\,\mathrm{d}\omega\\[6pt] &=\frac12\int_{-\infty}^\infty\frac1{\left(\omega_0^2-\omega^2\right)^2+(\Gamma_0\omega)^2}\,\mathrm{d}\omega\\ &=\frac1{4\omega_0^2}\int_{-\infty}^\infty\left[{\small\frac{1+\frac\omega{\sqrt{4\omega_0^2-\Gamma_0^2}}}{\left(\omega+\sqrt{\omega_0^2-\tfrac14\Gamma_0^2}\right)^2+\tfrac14\Gamma_0^2} +\frac{1-\frac\omega{\sqrt{4\omega_0^2-\Gamma_0^2}}}{\left(\omega-\sqrt{\omega_0^2-\tfrac14\Gamma_0^2}\right)^2+\tfrac14\Gamma_0^2}}\right]\,\mathrm{d}\omega\\ &=\frac1{4\omega_0^2}\int_{-\infty}^\infty\left[{\small\frac{\frac12+\frac\omega{\sqrt{4\omega_0^2-\Gamma_0^2}}}{\omega^2+\tfrac14\Gamma_0^2} +\frac{\frac12-\frac\omega{\sqrt{4\omega_0^2-\Gamma_0^2}}}{\omega^2+\tfrac14\Gamma_0^2}}\right]\,\mathrm{d}\omega\\[6pt] &=\frac1{4\omega_0^2}\int_{-\infty}^\infty\frac{\mathrm{d}\omega}{\omega^2+\tfrac14\Gamma_0^2}\\[12pt] &=\frac\pi{2\omega_0^2\Gamma_0} \end{align}This results in the following form for
$$ \langle x^2 \rangle = \sqrt{\dfrac{\pi}{2}} \dfrac{\tilde{f}}{4m^2\omega_0^2\Gamma_0} $$Assuming that the system is in thermal equilibrium and the equipartition theorem applies to this system.
From the equipartition theorem $$ \dfrac{1}{2}k_BT_{c.m.} = \dfrac{1}{2}m\omega_0^2 \langle x^2 \rangle $$ we get:
$$ \dfrac{1}{2}k_BT_{c.m.} = \dfrac{1}{2}m\omega_0^2 \sqrt{\dfrac{\pi}{2}} \dfrac{\tilde{f}}{4m^2\omega_0^2\Gamma_0} = \sqrt{\dfrac{\pi}{2}} \dfrac{\tilde{f}}{4m\Gamma_0} $$Plugging this into the previous result
$$ S_{xx}(\omega) = \sqrt{\dfrac{\pi}{2}} \dfrac{1}{m^2} \dfrac{\tilde{f}}{(\omega_0^2-\omega^2)^2 + (\Gamma_0\omega)^2} $$we get:
In [ ]: