I'm reading this paper Density of states for Gaussian unitary ensemble, Gaussian orthogonal ensemble, and interpolating ensembles through supersymmetric approach. In this work, the author developed a approach avoiding the Hubbard-Stratonovich transformation. To better understand my question, let me firstly introduce this approach in brief.
In terms of conventional and Grassmann numbers, we know the averaged resolvent can be written as $$ \frac{1}{N}\langle {\rm tr}\ G(E)\rangle=\frac{i}{N}\int\int[\boldsymbol{z},\boldsymbol{z}]\exp\left(-iE\left([\boldsymbol{z},\boldsymbol{z}]+[\boldsymbol{\psi},\boldsymbol{\psi}]\right)-\frac{1}{2N}\left([\boldsymbol{z},\boldsymbol{z}]^2-[\boldsymbol{\psi},\boldsymbol{\psi}]^2+2\sum_{jk}z_j\bar{\psi}_j\bar{z}_k\psi_k\right)\right) D\boldsymbol{z}D\boldsymbol{\psi}, $$ where $[\boldsymbol{z},\boldsymbol{z}]=\sum_j\bar{z}_jz_j=|\boldsymbol{z}|^2$ and the same fashion for Grassmann numbers $\boldsymbol{\psi}$, the covariance of the entries is set to $1/N$. The key point is the unitary transformation $U$ with $U_{1j}=\bar{z}_j/|\boldsymbol{z}|$. Denote $$ \bar{\phi}_i=\sum_j\bar{U}_{ij}\bar{\psi}_j,\ \phi_i=\sum_jU_{ij}\psi_j. $$ Then $$ [\boldsymbol{\phi},\boldsymbol{\phi}]=[\boldsymbol{\psi},\boldsymbol{\psi}],\ \sum z_j\bar{\psi}_j=|\boldsymbol{z}|\bar{\phi}_1,\ \sum \bar{z}_j\psi_j=|\boldsymbol{z}|\phi_1. $$ One can see that the cross term in the action $\sum_{jk}z_j\bar{\psi}_j\bar{z}_k\psi_k=[\boldsymbol{z},\boldsymbol{z}]\bar{\phi}_1\phi_1$ is decoupled after integrating over $\bar{\phi}_1,\phi_1$, which leads to $$ \frac{1}{N}\langle {\rm tr}\ G(E)\rangle=\frac{i}{N}\int\int[\boldsymbol{z},\boldsymbol{z}]\left(iE+\frac{1}{N}[\boldsymbol{z},\boldsymbol{z}]-\frac{1}{N}[\boldsymbol{\widetilde{\phi}},\boldsymbol{\widetilde{\phi}}]\right)\exp\left(-iE\left([\boldsymbol{z},\boldsymbol{z}]+[\boldsymbol{\widetilde{\phi}},\boldsymbol{\widetilde{\phi}}]\right)-\frac{1}{2N}\left([\boldsymbol{z},\boldsymbol{z}]^2-[\boldsymbol{\widetilde{\phi}},\boldsymbol{\widetilde{\phi}}]^2\right)\right) D\boldsymbol{z}D\boldsymbol{\widetilde{\phi}}, $$ where $\boldsymbol{\widetilde{\phi}}$ is the Grassmanian vector $\boldsymbol{\phi}$ without the first coordinate. Using the lemma proposed in the paper $$ \int D\boldsymbol{\phi}F\left([\boldsymbol{\phi},\boldsymbol{\phi}]\right)=(-1)^NF^{(N)}(0)=(-1)^N\frac{N!}{2\pi i}\oint \frac{F(z)}{z^{N+1}}dz, $$ where $F$ is any analytic function of the $N$-component Grassmann vector $\boldsymbol{\phi}$, then passing to the polar coordinates of $\boldsymbol{z}$ with the formula $$ \int_{\mathbb{C}^N}D\boldsymbol{z}F([\boldsymbol{z},\boldsymbol{z}])=\frac{1}{(N-1)!}\int_0^{\infty}F(s)s^{N-1}ds, $$ the final expression for the resolvent is $$ \frac{1}{N}\langle {\rm tr}\ G(E)\rangle=\frac{(-1)^{N-1}}{2\pi N}\int ds\oint dz\left(iE+s-z\right)e^{-N\left(f(s)+g(z)\right)}, $$ where $f(s)=s^2/2+iEs-\ln s$ and $g(z)=-z^2/2+iEz+\ln z$.
Question
In the large $N$ limit, I can finish the last integration with the sprit of steepest descent (saddle point) method. Now the confusing thing comes to me that if I directly put the saddle pair $(s_+,z_-)$ with $s_+=-iE/2+\sqrt{1-E^2/4}$ and $z_-=-s_+$ for $|E|<2$ into that integral, I will obtain $\Im\langle {\rm tr}\ G(E)\rangle/N<0$ for $E>0$, i.e the negative density of states.
Instead of the last integration, the author rewrited the resolvent as $$ \frac{1}{N}\langle {\rm tr}\ G(E)\rangle=is_++\frac{(-1)^{N-1}}{2\pi N}\int ds\oint dz\frac{iE+s-z}{s}(s-s_+)e^{-N\left(f(s)+g(z)\right)}. $$ This indeed covers the original expression for resolvent since the integral $$ \frac{(-1)^{N-1}}{2\pi N}\int ds\oint dz\frac{iE+s-z}{s}e^{-N\left(f(s)+g(z)\right)} $$ gives nothing but $i$ (two determinants cancels), thus gives the correct density of states. But how these two equivalent integration give different results, or in other word, what do I miss if I directly use the saddle point method in the last integral?
