$\mathbf{J}$ is a random matrix where $J_{ij}$ follows a Gaussian distribution.
Consider the following integral:
$$I=\int\left(\prod_{ij}\mathrm{d}J_{ij}\right) \exp\left\{-\frac{N}{2} \sum_{i, j, k} J_{k i} A_{i j} J_{k j}+N\sum_{k, j} B_{k j} J_{k j}\right\}$$
Where $\mathbf{A}$ and $\mathbf{B}$ are Hermitian. This is a regular Gaussian integral and by completing the square I can obtain (if not mistaken?):
$$I=(2 \pi)^{\frac{N^2}{2}}(\operatorname{det} \mathbf{A})^{-N / 2} \exp \left\{\sum_{i,j,k}^{n} \frac{1}{2} B_{ki}\left( A^{-1}\right)_{i j} B_{jk}\right\}$$
However if the elements $J_{ij}$ are correlated and my integral $I$ now becomes:
$$I=\int\left(\prod_{ij}\mathrm{d}J_{ij}\right) \exp\left\{-\frac{N}{2} \sum_{i, j, k} J_{k i} A_{i j} J_{k j}+N\sum_{k, j} B_{k j} J_{k j} +\tau N\sum_{ij}J_{ij}J_{ji}\right\}$$
with $-1<\tau<1$.
How can I deal with the $\sum_{ij}J_{ij}J_{ji}$ terms?
Any remark or advice is always appreciated. Thanks.
We assume that we are integrating over Hermitian matrices. Completing the square gives \begin{equation}\begin{aligned} &I=(2\pi)^{\frac{N^2}2}(\det A)^{-\frac N2}\int\prod_{i,j}dJ_{ij}\exp\Bigg(-\frac N2\text{tr}\bigg[((A-2\tau1_N)J-B)^\dagger(A-2\tau1_N)^{-1}((A-2\tau1_N)J-B) - \frac N2B(A-2\tau1_N)^{-1}B\bigg]\Bigg) \end{aligned}\end{equation} Doing a linear shift in $J$ by $(A-2\tau 1_N)^{-1}B$ gives us \begin{equation}\begin{aligned} &=(2\pi)^{\frac{N^2}2}(\det A)^{-\frac N2}\int\prod_{i,j}dJ_{ij}\exp\Bigg(-\frac N2\text{tr}\bigg[J^\dagger(A-2\tau1_N)J - \frac N2B(A-2\tau1_N)^{-1}B\bigg]\Bigg). \end{aligned}\end{equation}
So we need now to evaluate \begin{equation}\begin{aligned} Z=\int\prod_{i,j}dJ_{ij}\exp\Bigg(-\frac N2\text{tr}\bigg[J^\dagger AJ\bigg]\Bigg). \end{aligned}\end{equation} Since $A$ is Hermitian, there exists a unitary $U$ such that $A=UDU^\dagger$, for $D=\text{diag}(\lambda_1,\dots,\lambda_N)$ and we assume $\lambda_i\in\mathbb R_{>0}$. We do the change of variables $U^\dagger MU= J$. With this change of variables \begin{equation}\begin{aligned} \text{tr}(JAJ)&=\sum_i\lambda_iM_{ii}^2+\sum_{i\neq j}(\lambda_i+\lambda_j)\left((M_{ij}^{(r)})^2+(M_{ij}^{(im)})^2\right)\\ &=\sum_i\lambda_iM_{ii}^2+2\sum_{i<j}(\lambda_i+\lambda_j)\left((M_{ij}^{(r)})^2+(M_{ij}^{(im)})^2\right), \end{aligned}\end{equation} where $M_{ij}^{(r)}$ is the real part of $M_{ij}$ and $M_{ij}^{(im)}$ is the imaginary part.
Since $\det U$ has determinant 1, this means we can write \begin{equation}\begin{aligned} Z&=\int\prod_{i,j}dM_{ij}\exp\Bigg(-\frac N2\bigg[\sum_i\lambda_iM_{ii}^2+2\sum_{i<j}(\lambda_i+\lambda_j)\left((M_{ij}^{(r)})^2+(M_{ij}^{(im)})^2\right)\bigg]\Bigg)\\ &=\frac{(2\pi / N)^{N^2/2}}{\sqrt{\det A}\prod_{i<j}(\lambda_i+\lambda_j)}. \end{aligned}\end{equation}
Going back to $I$, if we assume that the eigenvalues of $A$ are all greater than $2\tau$, then the integral is convergent and we obtain \begin{equation} I=(2\pi/\sqrt{N})^{N^2}\exp\left(\frac N2\text{tr}\left[B(A-2\tau1_N)^{-1}B\right]\right)\frac{1}{\sqrt{\det(A-2\tau 1_N)}(\det A)^{N/2}\prod_{i<j}(\lambda_i+\lambda_j)}. \end{equation}
I'm going to guess that the normalisation constant for the integral is incorrect. If the normalisation constant was $$ C=\left(\frac{2\pi}{N^2/2}\right)^{-N^2}\sqrt{\det A}\prod_{i<j}(\lambda_i+\lambda_j), $$ then the integral would be 1 at $\tau=0$.
As you pointed out, this is not the problem you had. You had $J_{ij}$ where $J$ is real valued, and the term you added was Tr$(J^2)$. Now we can decompose $J=J^{(s)} + J^{(a)}$ where $J^{(s)}$ is symmetric and $J^{(a)}$ is antisymmetric. Then \begin{align} \text{Tr}(J^2)&=\text{Tr}((J^{(s)})^2 + (J^{(a)})^2+2J^{(s)}J^{(a)})\\ &=\text{Tr}((J^{(s)})^2 + (J^{(a)})^2)\\ &=\text{Tr}((J^{(s)})^TJ^{(s)} - (J^{(a)})^TJ^{(a)})\\ \text{Tr}(JAJ^T)&=\text{Tr}(J^{(s)}AJ^{(s)} - J^{(a)}AJ^{(a)})\\ &=\text{Tr}(J^{(s)}AJ^{(s)} + (J^{(a)})^TAJ^{(a)}). \end{align} This uses the fact that the trace of an antihermitian matrix is zero.
Now change variables to $K$ where $K^{(s)}=J^{(s)}$ and $K^{(a)}=iJ^{(a)}$. Then $K$ is Hermitian and this should reduce to the problem at the top of my answer, except for the issue now that the pure imaginary part multiply $A$ has a negative sign. This means that this integral will actually diverge. This can be saved if $\tau$ is bigger than the absolute value of the eigenvalues of $A$.