Consider the following Gaussian density over $\mathbb{R}^{2^n}$ $$p_n(\underline{x}):=\frac{\exp(-\frac{1}{2n}\langle C_n^{-1}\underline{x},\underline{x}\rangle)}{\sqrt{(2\pi n)^{2^n} \det C_n}}$$ where $C_n$ is a $2^n\times2^n$ real symetric positive defined matrix with entries bounded in $[-1,1]$.

Is there a strategy to compute the following integrals :

$$I_n := \int_{\mathbb{R}^{2^n}} \frac{\exp\left(\frac{\lambda}{n} x_1 x_2- b\ (x_1+x_2)\right)}{\left(\sum_{i=1}^{2^n}\exp(-b\ x_i)\right)^2}\ p_n(\underline{x})\ d\underline{x}$$

$$J_n := \int_{\mathbb{R}^{2^n}} \frac{\exp(\frac{\lambda}{n} x_1^2 - 2 b\ x_1)}{\left(\sum_{i=1}^{2^n}\exp(-b\ x_i)\right)^2}\ p_n(\underline{x})\ d\underline{x} $$

where $\lambda,b>0$ ?

In particular I'm interested in coumputing

$$\lim_{n\to\infty} 2^n(2^n-1)\ I_n + 2^n\ J_n$$

under some assuption of convergence of $C_n$.