Mean square displacement of a Langevin equation with inertia Consider a 1D motion of a particle
$$\ddot{x}(t)=-\gamma\dot{x}(t)+\eta(t)$$
where $\langle \eta \rangle=0$ and $\langle \eta(t)\eta(t') \rangle = \tilde{D} \delta(t-t')$.
How can I obtain analytically the long time asymptotic limit of the mean square displacement, i.e,
$$\langle x(t)-x(0))^2 \rangle$$
 A: This is a standard piece of bookwork, by which I mean that almost any textbook that discusses the Langevin equation will go on to derive an expression for the mean-square displacement. However, there are several different ways of formulating the derivation, some clearer than others. Therefore I'll just give a sketch of the route that I prefer.
The formal solution of the stochastic differential equation is, abbreviating $\dot{x}=v$:
$$
v(t) = \exp(-\gamma t) v(0) + \int_0^t dt' \exp[-\gamma (t-t')] \eta(t')
$$
Multiply by $v(0)$ and ensemble average,
the term involving the random force is zero ($\eta$ is uncorrelated with $v$ at earlier times),
and we get the well-known expression for the velocity autocorrelation function
$$
\langle v(t) v(0) \rangle = \langle v^2 \rangle \exp(-\gamma t).
$$
Strictly speaking it is $\exp(-\gamma |t|)$ but I assume $t>0$.
Now we write
$$
x(t) - x(0) = \int_0^t dt_1 v(t_1)
$$
square both sides and ensemble average
$$
\langle [x(t) - x(0)]^2\rangle = \int_0^t d t_1 \int_0^t d t_2 \langle v(t_1) v(t_2) \rangle
$$
At this point we note that the integrand is symmetric to the exchange $t_1\leftrightarrow t_2$.
The double integral range is illustrated here.
It is a square, which we divide into a dark grey and a light grey triangle.

For every point $(t_1,t_2)=(a,b)$ in the dark triangle (black dot) there is a mirror image point $(t_1,t_2)=(b,a)$ in the light triangle (white dot) which makes an identical contribution to the integral.
Therefore we may replace the integral by one over the dark triangle,
where $t_2$ ranges from $0$ to $t_1$ rather than $t$,
and double the result. 
Hence,
noting that $\langle v(t_1) v(t_2) \rangle = \langle v(t_1-t_2) v(0) \rangle=\langle v^2 \rangle\exp[-\gamma (t_1-t_2)]$:
\begin{align*}
\langle [x(t) - x(0)]^2\rangle &= 
2\langle v^2 \rangle \int_0^t d t_1 \int_0^{t_1} d t_2 \exp[-\gamma (t_1-t_2)]
\\
&=
2\langle v^2 \rangle \int_0^t d t_1 \exp(-\gamma t_1) \int_0^{t_1} d t_2 \exp(+\gamma t_2)
\end{align*}
The integral over $t_2$ yields $\gamma^{-1}[\exp(\gamma t_1)-1]$.
Multiplying this by $\exp(-\gamma t_1)$ and integrating over $t_1$
gives the result
$$
\langle [x(t) - x(0)]^2\rangle = \frac{2\langle v^2\rangle}{\gamma}
\left( t - \frac{1}{\gamma}[1-\exp(-\gamma t)]\right)
$$
The asymptotic limit is 
$$
\langle [x(t) - x(0)]^2\rangle \rightarrow \frac{2\langle v^2\rangle}{\gamma}t
= \frac{2 k_{\text{B}}T}{m\gamma} t = 2Dt
$$
where I have introduced the temperature, Boltzmann's constant, the mass and 
the diffusion coefficient $D$,
all of which can be related to your parameters $\gamma$ and $\tilde{D}$ through the fluctuation-dissipation theorem (your $\tilde{D}$ is not $D$!).
