We construct an $N\times N$ matrix $J$ whose elements are drawn from Gaussian distribution with zero mean and variance $\frac{1}{N}$. Since we want to have different variances for different columns, so that we express the elements as $J_{ij}\sigma_j$ such that the matrix satisfies $\mathbb{Var}[J_{ij}\sigma_j] = \frac{\sigma_j^2}{N}$.
Given a Green function (analogous electrostatic potential) $$\Phi(\omega) = -\frac{1}{N}\mathbb{E}_J\big[ \ln\det((I\omega^* - J^T)(I\omega - J)) \big]$$, where $\omega\in \mathbb{C}$ and $J$ is an $N\times N$ random matrix. Since $J$ is positive semi-definite, by adding diagonal matrix $\epsilon\delta_{ij}$, determinant can be represented by a Gaussian integral over complex variables, we have $$ \Phi(\omega) = \frac{1}{N}\ln\mathbb{E}_J\big[ \int\big( \prod_i\frac{d^2 z_i}{\pi}\big)\exp\big\{ -\epsilon\sum_i|z_i|^2 - \sum_{i, j, k} z_i^*(\omega^*\delta_{ik}-J_{ik}^T)(\omega\delta_{kj}-J_{kj})z_j \big\} \big] $$
Then the potential $\phi$ is determined by
$$ \exp(-N\Phi) = \int\prod_{i = 1}^N \frac{dz_i^* dz_i}{\pi}\exp(-N\hat{Q}) $$ with the average over the distribution of $J$ values given by $$ \exp(-N\hat{Q}) = \int\prod_{i,j = 1}^N \sqrt{\frac{N}{2\pi}}dJ_{ij}\exp(-NQ)$$ where $$ Q = \sum_i\big( \frac{|\omega|^2}{\sigma_i^2} + \epsilon_i \big)\frac{z_i^*z_i}{N} + \frac{1}{2}\sum_{i, j, k} J_{ki}A_{ij}J_{kj} - \sum_{k, j}B_{kj}J_{kj} $$ and $$ A_{ij} = \frac{z_i^*z_j}{N} + \frac{z_j^*z_i}{N} + \delta_{ij}\text{ and } B_{kj} = \frac{\omega^*z_k^*z_j}{\sigma_kN}+ \frac{\omega z_j^*z_k}{\sigma_kN} $$
Question:
It confuses me how this later expression expends the previous formula for $\Phi(\omega)$ ? Could someone help me to show some detailed calculations to explain ?