I have another way that yields a simple result that may be easier to deal with. I will use matrix notation, so dot products will be ##\mathbf{k \cdot x} = \mathbf{k}^T\mathbf{x}## where the T superscript represents transpose.
Let ##\mathbf{r}## be a 4 x 1 vector that contains ##\rho_1## and ##\rho_2## stacked on top of each other, and ##\mathbf{k}## be the 2x1 vector wave-vector. Then the integral is of the form
\begin{eqnarray}
I & = & \int_{-\infty}^{\infty} d^4 \mathbf{r} \, e^{i 2 \pi \mathbf{k}^T B \mathbf{r}}\, e^{-\mathbf{r}^T A \mathbf{r}}
\end{eqnarray}
Where ##B## is a 2x4 matrix and ##A## is a positive definite symmetric 4x4 matrix. For your problem I believe these have a very special form. If ##I_2## is a 2x2 identity matrix, then in block-matrix notation I think your problem can be represented
\begin{eqnarray}
B & = & \left( \begin{array} & I_2 & -I_2 \end{array} \right) \\
A & = & \left( \begin{array} & a I_2 & b I_2 \\ b I_2 & c I_2 \end{array} \right) \\
\end{eqnarray}
The results below do not rely on these forms, but they are probably helpful in doing the algebra to get the final answer in a form you want.
Anyway, now we use the fact that ##A## is positive definite to do an eigen-decomposition, ##A = P \Lambda P^T##, where ##\Lambda## is the diagonal matrix of eigenvalues and ##P## is an orthogonal matrix of eigenvectors. Then the integral can be written,
\begin{eqnarray}
I & = & \int_{-\infty}^{\infty} d^4 \mathbf{r} \, e^{i 2 \pi \mathbf{k}^T B \mathbf{r}}\, e^{-\mathbf{r}^T P \Lambda P^T \mathbf{r}}
\end{eqnarray}
Now we do a change of variable to ##\mathbf{y} = \Lambda^{1/2} P^T \mathbf{r}##. Note that the Jacobian is ##\left| \Lambda^{1/2} P^T\right|=\left| \Lambda\right|^{1/2} =\left| A \right|^{1/2}##, where ##| \cdot |## represents the determinant. So we have
\begin{eqnarray}
I & = & \left| A \right|^{-1/2}\int_{-\infty}^{\infty} d^4 \mathbf{y} \, e^{i 2 \pi \mathbf{k}^T B P \Lambda^{-1/2} \mathbf{y}}\, e^{-\mathbf{y}^T \mathbf{y}} \\
& = & \left| A \right|^{-1/2}\int_{-\infty}^{\infty} d^4 \mathbf{y} \, e^{\mathbf{s}^T\mathbf{y} - \mathbf{y}^T \mathbf{y}}
\end{eqnarray}
where ##\mathbf{s} = i 2 \pi \Lambda^{-1/2} P^T B^T \mathbf{k}##. Now we do the same trick as before to get rid the of the term linear in ##\mathbf{y}## (by letting ##\mathbf{z}=\mathbf{y}-\mathbf{s}/2##) and end up with,
\begin{eqnarray}
I & = & \left| A \right|^{-1/2}e^{\mathbf{s}^T\mathbf{s}/4} \int_{-\infty}^{\infty} d^4 \mathbf{z} \, e^{ - \mathbf{z}^T \mathbf{z}} \\
& = & \pi^2 \left| A \right|^{-1/2}e^{\mathbf{s}^T\mathbf{s}/4}
\end{eqnarray},
Now, ##\mathbf{s}^T\mathbf{s} = -4 \pi^2 \mathbf{k}^T B P \Lambda^{-1/2}\Lambda^{-1/2} P^T B^T \mathbf{k}=-4 \pi^2 \mathbf{k}^T B A^{-1} B^T \mathbf{k}##. So we finally have
\begin{eqnarray}
I = \pi^2 \left| A \right|^{-1/2}e^{-\pi^2 \mathbf{k}^T B A^{-1} B^T \mathbf{k}}.
\end{eqnarray}
The fact that ##A## and ##B## have such special forms will hopefully make this expression reducible. For example forumulas in
https://en.wikipedia.org/wiki/Block_matrix
can be used to find the inverse of ##A## quite easily, and I suspect the determinant of ##A## is also not so bad.
Good luck,
Jason
