I am trying to evaluate the definite integral below where $a, b, c$ are real:
$$ \int_{-1}^1 \int_{-1}^1 e^{-(ax^2 + bxy + by^2)} ~ \mathrm{d} x ~ \mathrm{d} y $$
I know that this can be done easily when integrating over $\mathbb{R}^2$, one can do a change of variables that maps the quadratic to another one where the $xy$ term vanishes. Since it is a linear change of variables, the limits of integration of the new integral are still from $-\infty$ to $+ \infty$, the integral is separable, and we are done.
If I attempt to do this on the above integral, I clearly end up with a region of integration that is no longer a rectangle aligned with the $x$ and $y$ axis, with the consequence the the inner integral's limits of integration depend on the outer variable of integration, and it no longer seems possible to solve it using the standard Gaussian integral approach.
It's pretty clear that this integral has no elementary closed form, but can it be expressed in terms of special functions (like the exponential integral and error function) or is there no way to proceed without a symmetry to exploit to separate the two variables?
It can be expressed using erf function $$ {\frac {\sqrt {\pi}}{2\sqrt{a}}\int_{0}^{1} {\rm e}^{{-by^2} +\frac{b^2y^2}{4a}} \left( {\rm erf} \left({\frac {by+2\,a}{2\sqrt{a}}}\right) -{\rm erf} \left({\frac {by}{2\sqrt{a}}}\right) \right) {\rm d}y } $$ but I am not sure it's very helpful...