This micro post is devoted to elementary integral formulas.
Integral representation of the determinant. If $\Sigma$ is a positive-definite symmetric $n\times n$ real matrix then the normalization of the density of the multivariate Gaussian law $\mathcal{N}(0,\Sigma)$ is
\[
Z:=\int_{\mathbb{R}^n}\mathrm{e}^{-\frac{1}{2}\langle\Sigma^{-1}x,x\rangle}\mathrm{d}x
=\sqrt{(2\pi)^n\det(\Sigma)}.
\] Conversely, this can also be used as an integral representation of the determinant or equivalently its logarithm : if $\Sigma$ is a positive-definite symmetric $n\times n$ real matrix then \[
\log\det(\Sigma)
=n\log(2\pi)
-2\log\int_{\mathbb{R}^n}\mathrm{e}^{-\frac{1}{2}\langle \Sigma x,x\rangle}\mathrm{d}x.
\]
Trace of the resolvent. If $S$ is an $n\times n$ symmetric matrix, and if $\lambda>-\min\mathrm{spec}(S)$, then the symmetric matrix $S+\lambda\mathrm{Id}$ is positive-definite and
\begin{align*}
\partial_\lambda\log\det(S+\lambda\mathrm{Id})
&=\sum_{i=1}^n\partial_\lambda\log(\lambda_i+\lambda)\\
&=\sum_{i=1}^n\frac{1}{\lambda_i+\lambda}\\
&=\mathrm{trace}((S+\lambda\mathrm{Id})^{-1}).
\end{align*}
Alternatively, with $V_\lambda(x):=\frac{1}{2}\langle(S+\lambda\mathrm{Id})x,x\rangle=\frac{1}{2}\langle\Sigma_\lambda^{-1}x,x\rangle$, and $X_\lambda\sim\mathcal{N}(0,\Sigma_\lambda)$,
\begin{align*}
-2\partial_\lambda\log\int_{\mathbb{R}^n}\mathrm{e}^{-V_\lambda(x)}\mathrm{d}x
&=-2\frac{\int_{\mathbb{R}^n}\partial_\lambda V_\lambda(x)\mathrm{e}^{-V_\lambda(x)}\mathrm{d}x}{\int_{\mathbb{R}^n}\mathrm{e}^{-V_\lambda(x)}\mathrm{d}x}\\
&=\frac{1}{Z_\lambda}\int_{\mathbb{R}^n}\|x\|^2\mathrm{e}^{-\frac{1}{2}\langle\Sigma_\lambda^{-1}x,x\rangle}\mathrm{d}x\\
&=\mathbb{E}[\|X_\lambda\|^2]\\
&=\mathrm{trace}(\Sigma_\lambda)\\
&=\mathrm{trace}((S+\lambda\mathrm{Id})^{-1}).
\end{align*} This gives also an integral representation of the trace of the resolvent \[
\mathrm{trace}((S+\lambda\mathrm{Id})^{-1})
=\frac{1}{Z_\lambda}
\int_{\mathbb{R}^n}\|x\|^2\mathrm{e}^{-\frac{\lambda}{2}\|x\|^2-\frac{1}{2}\langle Sx,x\rangle}\mathrm{d}x.
\] By a simple scale change of variable, we get, with $z=-\lambda$ and $\lambda$ large enough,\[
\mathrm{trace}((S-z\mathrm{Id})^{-1})
=\frac{1}{\widetilde{Z}_z}
\int_{\mathbb{R}^n}\|x\|^2\mathrm{e}^{-\frac{1}{2}\|x\|^2+\frac{1}{2z}\langle Sx,x\rangle}\mathrm{d}x.
\] Following Rǎzvan Gurǎu, beyond symmetric matrices, by replacing $\langle Sx,x\rangle$ by $\sum_{i_1,\ldots,i_p}T_{i_1,\ldots,i_p}$ for a $p$-fold tensor $T$, as soon as $z\in\mathrm{i}\mathbb{R}$, we get a trace of resolvent for symmetric tensors. This leads via series expansions to get a tensor analogue of the trace of powers. Alternatively, we could drop the restriction on $z$ and replace $\|x\|^2\mathrm{e}^{-\frac{1}{2}\|x\|^2}$ by $\|x\|^p\mathrm{e}^{-\frac{1}{2}\|x\|^p}$ to ensure integrability. Rémi Bonnin is working on this topic for his PhD.