Solution (source code)

= Solution

There is a source error to resolve first. For a <drag coefficient> defined by $F_{\mathrm{drag}}=-\xi v$, the <velocity relaxation time> is $\tau=m/\xi$. The mass and drag coefficient in the PDF's exponential are interchanged: $\xi t/m$ is dimensionless, whereas $mt/\xi$ is not. The consistent mean relaxation law is therefore $\langle v(t)\rangle=v_0e^{-t/\tau}$.

A mean relaxation law alone does not determine a <variance>. Use the usual thermal <Langevin equation>, with a bath at <temperature> $T$ and the <fluctuation-dissipation relation for a Langevin particle>:
$$
m\dot v=-\xi v+f(t),\qquad
\langle f(t)\rangle=0,\qquad
\langle f(t)f(t')\rangle=2\xi k_BT\,\delta(t-t').
$$
Here the calculation is for one component in <Cartesian coordinates>. With fixed initial position and velocity, integration gives the <mean displacement>
$$
\langle x(t)-x_0\rangle=v_0\tau(1-e^{-t/\tau}).
$$
Integrating the stochastic part of the velocity once more gives the centered displacement
$$
x(t)-\langle x(t)\rangle
=\frac1\xi\int_0^t\bigl[1-e^{-(t-s)/\tau}\bigr]f(s)\,ds.
$$
The delta correlation reduces its <variance> to a single <integral>:
$$
\begin{aligned}
\operatorname{Var}x(t)
&=\frac{2k_BT}{\xi}\int_0^t\bigl[1-e^{-(t-s)/\tau}\bigr]^2ds\\
&=\boxed{2D\left[t-2\tau(1-e^{-t/\tau})+\frac\tau2(1-e^{-2t/\tau})\right]},
\qquad \boxed{D=\frac{k_BT}{\xi}}.
\end{aligned}
$$
This is the <integrated Ornstein-Uhlenbeck displacement> for a prescribed initial velocity. Its short-time variance is $2Dt^3/(3\tau^2)+O(t^4)$; at long times it is $2Dt+O(1)$. Thus the <diffusion coefficient> is
$$
\boxed{D=\lim_{t\to\infty}\frac{\operatorname{Var}x(t)}{2t}}.
$$
If the requested mean-squared fluctuation instead refers to the conventional equilibrium <mean squared displacement>, the initial velocity is also averaged over its thermal distribution, with $\langle v_0\rangle=0$ and $\langle v_0^2\rangle=k_BT/m$. Its contribution adds $D\tau(1-e^{-t/\tau})^2$, giving
$$
\boxed{\langle[x(t)-x(0)]^2\rangle_{\mathrm{eq}}=2D\left[t-\tau(1-e^{-t/\tau})\right]}.
$$
This has the ballistic limit $(k_BT/m)t^2$ and the same diffusive limit $2Dt$. For fixed nonzero $v_0$, the uncentered mean square instead equals the first boxed <variance> plus $v_0^2\tau^2(1-e^{-t/\tau})^2$. In three independent spatial dimensions, the centered variances add, and the long-time result is $\langle|\mathbf r-\langle\mathbf r\rangle|^2\rangle\sim6Dt$. Stating the initial-velocity convention is necessary to distinguish these finite-time formulas.