Explicit solution for a linear drift-diffusion equation (Fokker-Planck equation) on whole space

I'm wondering if there might be an explicit solution for the following linear PDE in two space dimensions $$(x_1,x_2)$$ on the whole space $$\mathbb{R}^2$$: $$\partial_t f = {div} \left [\left( \begin{array}{rr} 1/4 & 0 \\ 0 & 1 \\ \end{array}\right) \nabla f + \left( \begin{array}{rr} 1/4 & -4 \\ 4 & 1 \\ \end{array}\right)xf \right ],$$ with $$x = (x_1,x_2) \in \mathbb{R}^2,$$ $$t >0$$ and $$f=f(t,x_1,x_2).$$ For this equation it is known that there is sufficient regular and positive, unique solution if we impose a positive initial state $$f_0 \in L^1(\mathbb{R}^2).$$ The problem is similar to the diffusion equation, for which, as is well known, there is a fundamental solution on the whole space. The equation can also be written as: $$\partial_t f = \frac{1}{4} \partial_{x_1 x_1}^2f + \partial_{x_2 x_2}^2f + \left(\frac{1}{4}x_1 - 4x_2 \right) \partial_{x_1} f + \left(4x_1 + x_2 \right) \partial_{x_2} f + \frac{5}{4}f.$$ The (normalized) steady-state $$f_{\infty}$$, solution of the corresponding elliptic equation $$0 = {div} \left [\left( \begin{array}{rr} 1/4 & 0 \\ 0 & 1 \\ \end{array}\right) \nabla f + \left( \begin{array}{rr} 1/4 & -4 \\ 4 & 1 \\ \end{array}\right)xf \right ]$$ is known, it would be $$f_{\infty}(x_1,x_2) = \frac{1}{2 \pi} {exp}\left(-\frac{1}{2}(x_1^2+x_2^2)\right).$$ Unfortunately I am not that very well versed in PDEs. So are there any results on explicit solutions of such equations on the whole space in 2D? I strongly assume the solution should be very similar to one of the heat equation.

I would be very grateful for any help!

The solution to the linear PDE is just the PDF corresponding to $$\mathcal{N}(\mu_t, \Sigma_t)$$ where $$\mu_t = \exp(A t) x \;, \quad \text{and} \quad \Sigma_t = 2 \int_0^t \exp(A (t-s) \begin{bmatrix} 1/4 & 0 \\ 0 & 1 \end{bmatrix} \exp(A^T (t-s)) ds$$ where $$A= -\begin{bmatrix} 1/4 & -4 \\ 4 & 1 \end{bmatrix}$$.
This solution can be easily obtained by using the correspondence between SDE/PDE. In particular, the corresponding SDE is $$d X_t = A X_t + \sqrt{2} \begin{bmatrix} 1/2 & 0 \\ 0 & 1 \end{bmatrix} d B_t$$ where $$B_t$$ is a standard $$2$$-dimensional Brownian motion. This SDE is a non-symmetric Ornstein-Uhlenbeck process due to the presence of a skew term in the drift, but since the SDE is linear, it can be explicitly solved to obtain the transition distribution $$\mathcal{N}(\mu_t, \Sigma_t)$$. It is also straightforward to verify that $$\mathcal{N}(\mu_t, \Sigma_t) \to \mathcal{N}(0,I_{2\times2})$$ as $$t \to \infty$$.
Metafune, G.; Pallara, D.; Priola, E., Spectrum of Ornstein-Uhlenbeck operators in $$L ^{p}$$ spaces with respect to invariant measures, J. Funct. Anal. 196, No. 1, 40-60 (2002). ZBL1027.47036.