Let $\mathrm{sn}(\theta) = \mathrm{sn}(\theta,\sqrt{m})$,
$\mathrm{cn}(\theta) = \mathrm{cn}(\theta,\sqrt{m})$,
$\mathrm{dn}(\theta) = \mathrm{dn}(\theta,\sqrt{m})$ be the
Jacobi elliptic functions associated with $k = \sqrt{m} = \sqrt{n(2-n)}$. Introduce parametrization:
$$
\theta = \int_{0}^t \frac{dt}{\sqrt{(1-t^2)(1-k^2t^2)}} \rightarrow t = \mathrm{sn}(\theta )
$$
and let $\mathrm{sn}(\theta_0)^2 = \frac{1}{n}$, we can rewrite the integral as:
$$
\Pi(n;\phi|k^2)= \mathrm{sn}(\theta_0)^2 \int_0^{\mathrm{sn}^{-1}(\sin\phi)} \frac{d\theta}{\mathrm{sn}(\theta_0)^2-\mathrm{sn}(\theta)^2}
$$
For given $\epsilon$, introduce $\varphi$, $\psi$ such that
$\sin(\varphi) = \sqrt{\frac{1-\epsilon}{n}} = \mathrm{sn}(\psi)$.
Notice
$$\begin{align}
\mathrm{cn}(\theta_0) &= \sqrt{1 - \mathrm{sn}(\theta_0)^2} = \sqrt{\frac{n-1}{n}}\\
\mathrm{dn}(\theta_0) &= \sqrt{1 - k^2\mathrm{sn}(\theta_0)^2} = \sqrt{n-1}\\
\frac{d}{d\theta} \mathrm{sn}(\theta) &= \mathrm{cn}(\theta)\mathrm{dn}(\theta)\\
\end{align}
$$
We see $\mathrm{cn}(\theta_0) = \mathrm{sn}(\theta_0)\mathrm{dn}(\theta_0)$ and
the expression appear within the definition of $C(n)$ becomes:
$$
\begin{align}
&2(n-1)\Pi\left(n;\varphi|k^2\right) +\log\epsilon\\
=& 2\,\mathrm{cn}(\theta_0)^2 \int_0^{\psi} \frac{d\theta}{\mathrm{sn}(\theta_0)^2-\mathrm{sn}(\theta)^2} + \log(1- \frac{\mathrm{sn}(\psi)^2}{\mathrm{sn}(\theta_0)^2} )\\
=& \int_0^{\psi} \frac{2 d\theta}{\mathrm{sn}(\theta_0)^2-\mathrm{sn}(\theta)^2}
( \mathrm{cn}(\theta_0)^2 - \mathrm{sn}(\theta)\mathrm{cn}(\theta)\mathrm{dn}(\theta) )\\
=& \int_0^{\psi} \frac{2 d\theta}{\mathrm{sn}(\theta_0)^2-\mathrm{sn}(\theta)^2}
( \mathrm{sn}(\theta_0)\mathrm{cn}(\theta_0)\mathrm{dn}(\theta_0) - \mathrm{sn}(\theta)\mathrm{cn}(\theta)\mathrm{dn}(\theta) )
\end{align}
$$
Since the numerator goes to $0$ as $\theta \to \theta_0$, the $\epsilon \to 0$ limit of the definition of $C(n)$ exists and is given by:
$$
C(n) = \int_0^{\theta_0} \frac{2 d\theta}{\mathrm{sn}(\theta_0)^2-\mathrm{sn}(\theta)^2}
( \mathrm{sn}(\theta_0)\mathrm{cn}(\theta_0)\mathrm{dn}(\theta_0) - \mathrm{sn}(\theta)\mathrm{cn}(\theta)\mathrm{dn}(\theta) )
$$
I don't know how to simplify this further but at least this is an alternate approach.
The reason you can't use the stationary phase technique here is that there's no point of stationary phase – the derivative of the phase
$$
\begin{align}
\phi(x)&=\delta x-\nu\exp(-\kappa x)/\kappa\;,
\\
\phi'(x)&=\delta+\nu\exp(-\kappa x)\;,
\end{align}
$$
is everywhere positive. You can turn that to your advantage though, since it allows you to change variables to $\phi$, yielding
$$
I=\int_{-\infty}^\infty \frac{G(x(\phi),\sigma)}{\phi'}\exp(-\mathrm i\phi)\,\mathrm d\phi\;.
$$
This is a Fourier component of a slowly varying function. For positive $x$ we have $\phi\sim\delta x$, whereas for negative $x$ we have $\phi\sim-\nu\exp(-\kappa x)/\kappa$ and thus $x/\sigma\sim\log(-\phi)/(\kappa\sigma)$, so the integral should decay exponentially with $\delta$ and $\kappa\sigma$. You can get upper bounds for its magnitude by integrating by parts; the boundary terms vanish because of the Gaussian, so you get
$$
|I|=\left|\int_{-\infty}^\infty \frac{G(x(\phi),\sigma)}{\phi'}\exp(-\mathrm i\phi)\,\mathrm d\phi\right|\le
\int_{-\infty}^\infty\left|\frac{\mathrm d^n}{\mathrm d\phi^n}\frac{G(x(\phi),\sigma)}{\phi'}\right|\,\mathrm d\phi
$$
for all $n\in\mathbb N$. Since each derivative yields a large inverse factor, this should lead to very small bounds on the magnitude of the integral.
These bounds are easier to evaluate if we transform back to $x$. For instance, for $n=0$ we have
$$
|I|\le\int_{-\infty}^\infty G(x,\sigma)\,\mathrm dx=1\;,
$$
for $n=1$, using $|\phi'|\lt\delta$ and $|\phi''/\phi'|\lt\kappa$,
$$
\begin{align}
|I|&\le\int_{-\infty}^\infty\left|\frac{G'(x(\phi),\sigma)}{\phi'^2}-\frac{G(x(\phi),\sigma)\phi''}{\phi'^3}\right|\,\mathrm d\phi
\\
&=\int_{-\infty}^\infty\left|\frac{G'(x,\sigma)}{\phi'}-\frac{G(x,\sigma)\phi''}{\phi'^2}\right|\,\mathrm dx
\\
&\le\frac1\delta\int_{-\infty}^\infty\left(|G'(x,\sigma)|+\kappa|G(x,\sigma)|\right)\,\mathrm dx
\\
&=\frac1\delta\left(2G(0,\sigma)+\kappa\right)
\\
&=\frac1\delta\left(\frac2{\sqrt{2\pi}\sigma}+\kappa\right)\;,
\end{align}
$$
and for $n=2$, using $|\phi'''/\phi'|\le\kappa^2$,
$$
\begin{align}
|I|&\le\int_{-\infty}^\infty\left|\frac{G''}{\phi'^3}-3\frac{G'\phi''}{\phi'^4}+3\frac{G\phi''^2}{\phi'^5}-\frac{G\phi'''}{\phi'^4}\right|\,\mathrm d\phi
\\
&=\int_{-\infty}^\infty\left|\frac{G''}{\phi'^2}-3\frac{G'\phi''}{\phi'^3}+3\frac{G\phi''^2}{\phi'^4}-\frac{G\phi'''}{\phi'^3}\right|\,\mathrm dx
\\
&\le\frac1{\delta^2}\int_{-\infty}^\infty\left(|G''|+3\kappa |G'|+3\kappa^2|G|+\kappa^2|G|\right)\,\mathrm dx
\\
&\le\frac1{\delta^2}\left(\frac2{\sigma^2}+\frac{6\kappa}{\sqrt{2\pi}\sigma}+4\kappa^2\right)\;.
\end{align}
$$
We can go on like this until the cows come home, that is, until the numerical factors from the derivatives begin to overwhelm the factors of $\kappa/\delta$, so if $\kappa\ll\delta$, the integral is for all intents and purposes zero. If $\kappa\sim\delta$, you may have to work a bit harder to get better bounds.
Best Answer
This is just an outline of an answer.
Call the integral to be calculated $I(a)$. First write $$I(a) = 2\int_0^{\pi/2} \sqrt{\cos^2 k + a^2 \sin^2 k} \, dk$$ by symmetry. The difficulty here is that you can't just apply Taylor's formula for $a \to 0$, because values of $k$ near $\pi/2$ make a significant contribution to the integral and $\cos k$ is small there.
Make the substitution $u = \cot k$. Then $$I(a) = 2 \int_0^{+\infty} \frac{\sqrt{u^2 + a^2}}{(u^2 + 1)^{3/2}} \, du.$$
Now divide the integral into three parts on the intervals $[0,a]$, $[a,1]$ and $[1,+\infty)$. Then make the substitution $v = u/a$ on the first interval, $s = u^2$ on the second, and $w = 1/u$ on the third. We get $$I(a) = 2a^2 \int_0^1 \sqrt{1+v^2} (1 + a^2 v^2)^{-3/2} \, dv + 2\int_0^1 (1 + w^2)^{-3/2} \sqrt{1 + a^2 w^2} \, dw \\+ \int_{a^2}^1 \frac{1}{(1+s)^{3/2}} \sqrt{1 + \frac{a^2}{s}} \, ds$$ Now the idea is to expand each integrand into a series using the Taylor series for $(1 + x)^{1/2}$ and $(1+x)^{-3/2}$, and then integrate term by term. This will be legitimate because of uniform convergence. In the first integral, expand the second factor as a function of $av$. Do the same in the second integral with respect to $aw$. The third is more complicated because $a^2$ appears as a bound, but you can expand the integrand into a double series with respect to $s$ and $a^2/s$. The resulting series is quite complicated.
However, if we only want an estimate at the level of $O(a^2)$, we can note that the square root in the last integral is $1 + a^2/2s$ to within $O(a^4/s^2)$, so the error in the integral will be at most $O(a^2)$. If we make this approximation, we find $$I(a) = 2 - a^2 \ln a + O(a^2)$$
To evaluate the $a^2$ term is more difficult. The contribution from the first integral is $\sqrt{2} + \operatorname{arsinh}(1)$. The contribution from the second is $\operatorname{arsinh}(1) - \frac{1}{2}\sqrt{2}$. The $a^2$ term in the third integral consists of a $-a^2$ from the first term of the series, a term $a^2[-\operatorname{arsinh}(1)+ \frac{1}{\sqrt{2}} + \ln 2 - 1]$ in the second term of the series, and the remaining terms with coefficient $\sum_{n \geq 2} \frac{1}{n-1}\binom{1/2}{n}$. This last series is $f(1)$, where $f(x) = \sum_{n \geq 2} \frac{1}{n-1}\binom{1/2}{n}x^{n-1}$. We have $f(0) = 0$ and $f'(x) = \sum_{n \geq 2} \binom{1/2}{n}x^{n-2} = \frac{1}{x^2} (\sqrt{1 + x} - 1 - x/2)$, so $$f(1) = \int_0^1 \frac{1}{x^2} (\sqrt{1 + x} - 1 - x/2) \, dx = \frac{3}{2} - \sqrt{2} + \ln 2 - \operatorname{arsinh}(1).$$ Taking everything into account, we get $$I(a) = 2 - a^2\ln a + a^2 (2\ln 2 - 1/2) + O(a^4 \ln a).$$
Given how simple the result is, I'll bet there's a simpler way to find it.
EDIT: For $a=0.0001$, the true value of the integral is $2.000,000,100,966,347,688$. The approximation obtained by the formula is $2.000,000,100,966,347,331$.