\begin{align} \int\int{\exp(-k|x_1-x_2|^2)d^2x_1d^2x_2} \end{align}
How to perform the integration of the integral above by changing the variable to the center of mass and relative coordinate,as below, \begin{align} u&=(x_2+x_1)\\ v&=\frac{1}{2}(x_2-x_1) \end{align}
On
I'll set $x_1=x$ and $x_2=y,$ and also use the familiar notation for differentials. Then I'll call your integrand $F(x,y)$ and your equations $x=f(u,v),\,y=g(u,v).$
Calculate the Jacobian of the transformation $(u,v)\mapsto (x,y).$ This is the determinant formed by the partials of $x=f(u,v),\,y=g(u,v).$ In particular, the first row contains the partials of $x$ wrt $u,v$ in that order, and the second row something similar for $y.$ Then your integral is $$\iint F(x,y)|J(u,v)|\mathrm du \mathrm dv,$$ where $x=f(u,v),\,y=g(u,v)$ from your equations, provided $J(u,v)$ does not vanish.
By the way you didn't indicate the region of integration, so this is as far as can be said.
This is about the quantity $$Q:=\int_{{\mathbb R}^2\times{\mathbb R}^2}\exp\bigl(-k|x-y|^2\bigr)\>{\rm d}(x_1,x_2,y_1,y_2)\ .$$ Unfortunately $Q=\infty$. In order to show this we use the transformation $$x_1={u_1+v_1\over2},\quad y_1={v_1-u_1\over2},\quad x_2={u_2+v_2\over2},\quad y_2={v_2-u_2\over2}\ .\tag{1}$$ Then $x-y=(u_1,u_2)$, and the Jacobian of $(1)$ computes to ${1\over4}$. It follows that $$Q={1\over4}\int_{{\mathbb R}^2\times{\mathbb R}^2}\exp\bigl(-k(u_1^2+u_2^2)\bigr)\>{\rm d}(u_1,u_2,v_1,v_2)\ .$$ Here the $u$-integral is fine, but the $v$-integral is $=\infty$.