Is there, similar to the Mehler kernel, a closed formula for the heat kernel of the heat equation associated to the Laplacian $$ \sum_j \frac{d^2}{dx_j^2} + 2\sqrt{1} \sum_j \lambda_j \frac{d}{dx_j} + \sum_{ij} a_{ij}x_ix_j$$ on $\mathbb{R}^n$? Here, the matrix $(a_{ij})$ is supposed to be symmetric and positive definite, while the $\lambda_j$ can be arbitrary.

Yes there is. Here is how you do it. First find an orthogonal change in variables $$ x_j=\sum_{jk} s_{jk}y_k $$ $(s_{jk})$ orthogonal matrix, so that in the new coordinates we have $$ \sum_{i,j}a_{ij} x_ix_j = \sum_j \mu_j^2 y_j^2, $$ where $\mu_j^2$ are the eigenvalues of the symmetric matrix $(a_{ij})$. Note that $\newcommand{\ii}{\boldsymbol{i}}$ $\newcommand{\pa}{\partial}$ $$\frac{\pa}{\pa x_k}=\sum_j\frac{\pa y_j}{\pa x_k}\frac{\pa }{\pa y_j} =\sum_j t_{jk}\frac{\pa}{\pa y_j}, $$ where $(t_{jk})$ is the inverse of the orthogonal matrix $(s_{jk})$ so that $t_{jk}=s_{kj}$. Then for some real numbers $r_j$ $$ \sum_j\frac{\pa^2}{\pa x_j^2}+2\ii\sum_j\lambda_j\frac{\pa }{\pa x_j}+\sum_{i,j}a_{ij}x_ix_j $$ $$ =\sum_j\frac{\pa^2}{\pa y_j^2} +2\ii\sum_j r_j\frac{\pa}{\pa y_j} +\sum_j \mu_j^2 y_j^2 $$ $$=\underbrace{\sum_j \left(\ii\frac{\pa}{\pa y_j}+r_j\right)^2 +\sum_j\mu_j^2y_j^2\sum_j r_j^2}_{=: L}. $$ Next set $$ R^2 :=\sum_j r_j^2, \;\; w(t,y) := R^2t +\sum_j \ii r_j y_j, $$ $$ L_0 :=\sum_j\left(\frac{\pa^2}{\pa y_j^2} +\mu_j^2y_j^2\right), $$ and observe that $$\pa_t +L =e^{w}(\pa_t+L_0) e^{w}. $$ Thus, if $K$ is a fundamental solution of $\pa_t+L_0$, then $$(\pa_t +L) (e^{w} K) = e^{w} (\pa_t+L_0)K=e^{w} \delta_0=\delta_0 $$ so that $e^{w}K$ is a fundamental solution of $\pa_t+L$. You can find $K$ using Mehler's formula. 


Lars Hörmander did some work on classifying Mehlertype formulas for general quadratic fomrs in $\xi$ and $x$. Take a look at his paper on Math.Zeitschrift 219 (1995) pp.413449 

