Take the stochastic process coming from
$$\ddot X + \frac{\eta}{m} \dot X + \frac{k}{m}X = \frac{1}{m} W(t)\tag{1}$$
(where $W$ is white noise) or equivalently the Itô differential equation
$$ d \begin{pmatrix} X\\V\\\end{pmatrix} = \begin{pmatrix} V\\ -\dfrac{\eta}{m} V - \dfrac{k}{m}X\\\end{pmatrix}dt + \begin{pmatrix} 0&0\\0&\frac1m\end{pmatrix}\begin{pmatrix} d\beta_1 \\d\beta_2\end{pmatrix}\tag{2} $$
$$ \equiv d \mathbf{X} = f(\mathbf{X})dt + \mathbf{L}\, d\boldsymbol{\beta} \tag{3}.$$
This is just the damped harmonic oscillator driven by white noise. In particular, $\beta_1,\beta_2$ are Brownian motion. I know from the literature that the stationary distribution is of the form
$$p(x,v) = \exp( - Ax^2 - Bv^2 )$$ (see edit).
Attempt: I believe I need to use the Fokker-Planck equation:
$$\frac{\partial p(\mathbf{X},t)}{\partial t} = 0= - \sum_{i}\frac{\partial}{\partial X_i}[\mathbf{f}_i(\mathbf{X})p(\mathbf{X},t)]+ \sum_{i,j} \frac{\partial^2}{\partial X_i X_j}\left[(\mathbf{L}^T\mathbf{L})_{ij} \,p(\mathbf{X},t)\right]\tag{4}$$ where the $0$ is because we consider the stationary distribution. This, however, gives me (writing $p(\mathbf{X},t) = p(x,v)$)
$$0=-\left[ \frac{\partial}{\partial x}(vp(x,v)) + \frac{\partial}{\partial v}\left(\left( -\frac{\eta}{m}v - \frac{k}{m}x\right) p(x,v)\right) \right]+ \frac12 \frac{\partial^2}{\partial v^2} \frac{p(x,v)}{m^2}\tag{5}$$
which is certainly not solved by a Gaussian $p(x,v)$.
Edit: An error stopped me from solving (5) with a Gaussian ansatz and plugging and chugging. Thank you stochastic for clearing this up.
I have confirmed that $(2)$ has the following distribution: $$P(x,v)=\frac{\eta\sqrt{km}}{\pi}e^{-\eta(kx^2+mv^2)}.$$
It also satisfies $(5).$
Details: Let us rewrite your equation in a matrix form: $$d\vec x= A\vec x\, dt+L\,d\vec \beta,$$ where $$A = \left(\begin{array}{cc}0&1\\-\frac km&-\frac\eta m\end{array}\right).$$
First, notice that starting with a Gaussian distribution, the right-hand side of the equation is a linear combination of Gaussians, and there the change in $x$ will be Gaussian and therefore $x$ stays Gaussian. So the steady-state solution should be a Gaussian. All you need to do is find the covariance $\newcommand\mean[1]{\left\langle#1\right\rangle} C=\mean{\vec x .\vec x^\top}.$ Using Ito's lemma we have: $$\frac{dC}{dt}=\frac{d}{dt}\mean{\vec x .\vec x^\top}=A\mean{\vec x .\vec x^\top}+\mean{\vec x .\vec x^\top}A^\top+L^\top L = AC+CA^\top+L^\top L. $$
At steady state, $AC+CA^\top+L^\top L=0$ can be solved for $C$: $$C=\left(\begin{array}{cc}\frac1{2k\eta}&0 \\ 0&\frac1{2\eta m}\end{array}\right).$$ Given the covariance matrix the Gaussian distribution is given by $$ P(\vec x) = \frac1{\sqrt{det(2\pi C)}}e^{-\frac12\vec x^\top.C^{-1}.\vec x}. $$