A formula is given by Wiki (https://en.wikipedia.org/wiki/Gaussian_integral#n-dimensional_with_linear_term):
$$ \int_{\mathbb{R}^n} \mathrm{d}^n\mathbf{x}\ \exp(-\frac{1}{2}\mathbf{x}^T\mathbf{A}\mathbf{x}+\mathbf{j}^T\mathbf{x}) = \sqrt{\frac{(2\pi)^n}{\det\mathbf{A}}}\exp(\frac{1}{2}\mathbf{j}^T\mathbf{A}^{-1}\mathbf{j}) $$ where $\mathbf{j}$ is a vector. It remains unclear to me what are the constraints on $\mathbf{j},\mathbf{A}$? Clearly, if both are real-valued and $\mathbf{A}$ is positive definite the the identity holds. Does it hold for complex $\mathbf{j}$ and Hermitian positive definite $\mathbf{A}$ (while integrating over $\mathbb{R}^n$ only)?
A few numerical tests suggest so and I can't see why it would not.