Reference for permanent integral identity $\DeclareMathOperator\perm{perm}\DeclareMathOperator\diag{diag}$Using MacMahon's master theorem, the properties of complex gaussian integrals, and Cauchy's integral theorem one can show that the permanent of a matrix $A$ satisfies
\begin{equation}
\perm(A) = \frac{1}{\pi^N} \int_{\mathbb{R}^{2N}} d^N\mathbf{u} \, \exp\Bigl( -\sum_{i=1}^N \mathbf{u}_i^2 \Bigr) \prod_{i=1}^{N} \Bigl(\mathbf{u}_i\cdot\sum_{j=1}^NA_{ij} \mathbf{u}_j\Bigr),
\end{equation}
where $\mathbf{u}_i \equiv (u_{i1}, u_{i2})$ and $d^{N}\mathbf{u} \equiv \prod^{N}_{i =1}du_{i1} \wedge du_{i2}$.
My guess is that this result is not new, but I'm not familiar with the literature on such identities. Is this identity known? Alternatively, where would be a good place to check for similar such identities?

Proof of Identity
The truth of the identity could be inferred from noting that we need our gaussian to integrate a function where $\mathbf{u}_i$ for each $i$ appears exactly twice.
To prove it directly, we note that MacMahon's master theorem states that the permanent of a matrix $A$ is the $x_1x_2\dotsm x_N$ coefficient of the quantity
\begin{equation}
\frac{1}{\det(I - XA)},
\end{equation}
where $I$ is the $N\times N$ identity matrix and $X \equiv \diag(x_1, x_2, \dotsc, x_N)$. Thus, we have
\begin{equation}
\perm(A) = \frac{1}{(2\pi i)^N} \oint \left[\prod_{i=1}^{N} \frac{dq_i}{q_i^2} \right] \frac{1}{\det(I- QA)}, 
\label{eq:perm_def}
\end{equation}
where $Q = \diag(q_1, q_2, \ldots, q_N)$, $q_i$ is a complex variable, and we are performing $N$ contour integrations in sequence. Defining a set of $N$ two-dimensional vectors $\{\mathbf{u}_i\}$ as $\textbf{u}_i \equiv (u_{i1}, u_{i2})$, we have
\begin{equation}
\frac{1}{\det(I-QA)} = \frac{1}{\pi^N} \int_{\mathbb{R}^{2N}} d^{N}\mathbf{u}\,\exp\Bigl(- \sum_{i, j=1}^N \mathbf{u}_i \cdot \left(\delta_{ij} - q_{i} A_{ij}\right) \mathbf{u}_{j}\Bigr).
\label{eq:gauss_result}
\end{equation}
Using this determinant expression in MacMahon's theorem and performing the contour integrations yields the stated result.
 A: In case anyone from the future comes to this post:
I wasn't able to find relevant literature that answered my question so I wrote a preprint on it (https://arxiv.org/abs/2106.11861). It turns out that the gaussian above is more complicated than necessary. The main result of the paper is a generalization of the identity above:

Let $p_X: \Omega_X \to \mathbb{R}$ be a probability distribution defined over the domain $\Omega_X$ with zero mean and unit variance. Let $A$ be an $n\times n$ matrix with elements $a_{i, j}$. Then the permanent of $A$ is
\begin{equation}
\text{perm}(A) = \int_{\Omega^n_{X}} d^n\textbf{x}\, \prod_{i=1}^n p_{X}(x_i) \, x_{i} \sum_{j=1}^n a_{i, j} x_j,
\label{eq:fund_thm}
\end{equation}
where $\Omega^n_X = \Omega_X  \otimes \cdots \otimes \Omega_X$ is the $n$-factor product over the single-variable domain of integration. In condensed notation, we can write this result as the expectation value
\begin{equation}
\text{perm}(A) = \left\langle \prod_{i=1}^nx_{i} \sum_{j=1}^n a_{i, j} x_j \right\rangle_{x_i \sim p_X},
\end{equation}
where the average is over $\{x_i\}$, a set of independent identically distributed random variables each of which is drawn from $p_X$.

When $p_X = N(0, 1)$, we get something like the equation in the question.
