Integral $\iint xy\frac{1}{\sqrt{b}}\exp\left[-\frac{(x-y)^2}{2b}\right]\sqrt{\frac{2}{c}}\exp\left[-\frac{(x+y)^2}{4c}\right].$ I am interested in the following integral 
$$\int_{-a}^adx\int_{-a}^a\mathop{\mathrm{d}y}xy\frac{1}{\sqrt{b}}\exp\left[-\frac{(x-y)^2}{2b}\right]\sqrt{\frac{2}{c}}\exp\left[-\frac{(x+y)^2}{4c}\right].$$
Does any one know how to evaluate the above integral? Mathematica was not able to give result for this. 
 A: Let's make this integral a little nicer. First, we deal with the parameters:
$$\int_{-a}^a\int_{-a}^axy\frac{1}{\sqrt{b}}\exp\left[-\frac{(x-y)^2}{2b}\right]\sqrt{\frac{2}{c}}\exp\left[-\frac{(x+y)^2}{4c}\right]\mathop{\mathrm{d}x}\mathop{\mathrm{d}y}=$$

$$x=au, \qquad y=av, \qquad b=a^2 B, \qquad c=\frac{a^2 C}{2}$$

$$=\frac{2a^2}{\sqrt{BC}} \int_{-1}^1 \int_{-1}^1 u~v~\exp \left[-\frac{(u-v)^2}{2B} \right]~\exp \left[-\frac{(u+v)^2}{2C} \right]\mathrm{d}u~\mathrm{d}v$$

Let's forget about the coefficients and name the integral:
$$I_0(B,C)=\int_{-1}^1 \int_{-1}^1 u~v~\exp \left[-\frac{(u-v)^2}{2B} \right]~\exp \left[-\frac{(u+v)^2}{2C} \right]\mathrm{d}u~\mathrm{d}v$$
Separating the region into four obvious parts and introducing another integral:
$$I(B,C)=\int_0^1 \int_0^1 u~v~\exp \left[-\frac{(u-v)^2}{2B} \right]~\exp \left[-\frac{(u+v)^2}{2C} \right]\mathrm{d}u~\mathrm{d}v$$
We can easily see that:

$$I_0(B,C)=2I(B,C)-2I(C,B)$$

Thus, we only need to find $I(B,C)$. Let's expand the squares in the exponentials:
$$I(B,C)=\int_0^1 \int_0^1 u~v~\exp (-\alpha ~u~v )~\exp \left[-\beta (u^2+v^2) \right]\mathrm{d}u~\mathrm{d}v$$
$$\alpha=\frac{1}{C}-\frac{1}{B}, \qquad \beta=\frac{1}{2C}+\frac{1}{2B}$$
Let's introduce a new integral:
$$J(\alpha,\beta)=\int_0^1 \int_0^1 \exp (-\alpha ~u~v )~\exp \left[-\beta (u^2+v^2) \right]\mathrm{d}u~\mathrm{d}v$$

$$I(B,C)=-\frac{\partial J(\alpha,\beta)}{\partial \alpha}$$


Now we notice that the function under the integral is symmetric with respect to the change $u \to v$, thus we can write:
$$J(\alpha,\beta)=2\int_0^1 \int_0^u \exp (-\alpha ~u~v )~\exp \left[-\beta (u^2+v^2) \right]\mathrm{d}u~\mathrm{d}v$$
Now we make a change of variable:
$$v=u t, \qquad dv=u~dt$$
$$J(\alpha,\beta)=2\int_0^1 \int_0^1 \exp (-\alpha ~u^2~t )~\exp \left[-\beta u^2(1+t^2) \right]u~\mathrm{d}u~\mathrm{d}t$$
Another change of variable will be:
$$u^2=p$$
$$J(\alpha,\beta)=\int_0^1 \int_0^1 \exp [-(\beta t^2+\alpha ~t+\beta)p ]\mathrm{d}p~\mathrm{d}t$$
The inner integral is elementary and we obtain a single integral:


$$J(\alpha,\beta)=\int_0^1 \frac{1-\exp [-(\beta t^2+\alpha ~t+\beta) ]}{\beta t^2+\alpha ~t+\beta}~\mathrm{d}t$$
$$=\frac{1}{\beta} \int_0^1 \frac{\mathrm{d}t}{(t+\gamma)^2 +1-\gamma^2}-e^{\beta (1-\gamma^2)} \frac{1}{\beta} \int_0^1 \frac{\exp [-\beta (t+\gamma)^2 ]}{(t+\gamma)^2 +1-\gamma^2}~\mathrm{d}t$$
$$\gamma=\frac{\alpha}{2 \beta}=\frac{B-C}{B+C}$$


The first part is an elementary integral (arctangent if $bc>0$).
The second part is related to the error function.

Let's get back to the initial integral and express it in the new form. Note that the exchange $B \to C$ is equivalent to $\alpha \to - \alpha$ and $\beta \to \beta$:
$$\int_{-a}^a\int_{-a}^axy\frac{1}{\sqrt{b}}\exp\left[-\frac{(x-y)^2}{2b}\right]\sqrt{\frac{2}{c}}\exp\left[-\frac{(x+y)^2}{4c}\right]\mathop{\mathrm{d}x}\mathop{\mathrm{d}y}=$$

$$=2a^4\sqrt{\frac{2}{bc}} \left(-\frac{\partial J(\alpha,\beta)}{\partial \alpha}+\frac{\partial J(-\alpha,\beta)}{\partial \alpha} \right) =2a^4\sqrt{\frac{2}{bc}} \times$$
$$ \times \int_0^1 \left[ \frac{ 1-\left(\beta t^2 +\alpha  t+\beta+1\right) e^{-(\beta t^2 +\alpha  t+\beta)}}{(\beta t^2 +\alpha  t+\beta)^2}-\frac{ 1-\left(\beta t^2 -\alpha  t+\beta+1\right) e^{-(\beta t^2 -\alpha  t+\beta)}}{(\beta t^2 -\alpha  t+\beta)^2}~\right] t \mathrm{d}t $$

Numerically the first and last expressions are equal, which was checked with Mathematica for random values of $a,b,c$.
