Find the value of $\int_{-\infty}^\infty \int_{-\infty}^\infty e^{-(x^2+xy+y^2)} \, dx\,dy$

For fun, I want to point out that much more can be said. The spectral theorem for real symmetric matrices tells us that real symmetric matrices are orthogonally diagonalizable. Thus, if $A$ is some symmetric matrix, then there exists an orthogonal matrix $U$ and diagonal matrix $D$ such that $A$ can be written as $A=U^{-1} DU$. Hence $x^TAx=x^TU^TDUx=(Ux)^TD(Ux)$. The Jacobian det of the transformation $x\mapsto Ux$ is simply $1$, so we have a change of variables $u=Ux$:

$$\begin{align} \int_{{\bf R}^n}\exp\left(-x^TAx\right)\,dV & =\int_{{\bf R}^n}\exp\left(-(\lambda_1u_1^2+\cdots+\lambda_nu_n^2)\right)\,dV \\[6pt] & =\prod_{i=1}^n\int_{-\infty}^{+\infty}\exp\left(-\lambda_i u_i^2\right)du_i \\[6pt] & =\prod_{i=1}^n\left[\frac{1}{\sqrt{\lambda_i}}\int_{-\infty}^{+\infty}e^{-u^2} \, du\right] \\[6pt] & =\sqrt{\frac{\pi^n}{\det A}}.\end{align}$$

Note that we don't even have to compute $U$ or $D$. The above calculation is the generalized version of the "completing the squares" approach when $n=2$ (which is invoked elsewhere in this thread).

This formula is in fact the basis for the Feynman path integral formulation of the functional determinant from quantum field theory; since technically the integral diverges we need to compare them instead of looking at individual ones outright. Wikipedia has some more details.

$$ \begin{align} \int_{-\infty}^\infty\int_{-\infty}^\infty e^{-(x^2+xy+y^2)}\,\mathrm{d}x\,\mathrm{d}y &=\int_{-\infty}^\infty\int_{-\infty}^\infty e^{-((x+y/2)^2+3y^2/4)}\,\mathrm{d}x\,\mathrm{d}y\\ &=\int_{-\infty}^\infty\int_{-\infty}^\infty e^{-(u^2+3y^2/4)}\,\mathrm{d}u\,\mathrm{d}y\\ &=\int_{-\infty}^\infty e^{-u^2}\,\mathrm{d}u \int_{-\infty}^\infty e^{-3y^2/4}\,\mathrm{d}y\\ &=\int_{-\infty}^\infty e^{-u^2}\,\mathrm{d}u \sqrt{\frac43}\int_{-\infty}^\infty e^{-v^2}\,\mathrm{d}v\\ &=\sqrt\pi\sqrt{\frac43}\sqrt\pi\\ &=\pi\sqrt{\frac43}\\ &=\frac{2\pi}{\sqrt3} \end{align} $$