QHO in Microcanonical Ensemble: Problem with alternate derivation I am working through Franz Schwabl's book on Statistical Mechanics, and he has a number of derivations of thermodynamic quantities that are different than those I have seen before. I am also having difficulty finding them repeated elsewhere.
In particular, he has a method for calculating $\Omega(E)$, the number of states with a given energy $E$, of a series of $N$ independent Quantum Harmonic Oscillators ($\mathcal{H} = \sum_{j=1}^N\hbar\omega(n_j+\frac{1}{2})$) that I hadn't seen before. Proceeding from the result
$$\Omega(E) = \mathrm{Tr}\,\delta(\mathcal{H}-E)=\sum_{n_1=0}^{\infty}\cdots\sum_{n_N=0}^{\infty}\delta\left(E - \hbar\omega\sum_{j=1}^N\left(n_j+\frac{1}{2}\right)\right),$$
my strategy would be combinatoric: the delta-function turns the unrestricted sums over $n_j$ to a constraint on the total number of quanta. Calculating the number of ways you can partition $n=\sum_{j=1}^Nn_j$ quanta of energy among $N$ oscillators gives you $\Omega(E)$. This is the way we did it in undergraduate stat mech.
Schwabl's approach proceeds differently: by taking the Fourier Transform of the delta function, one obtains
$$\Omega(E) = \int \frac{dk}{2\pi}e^{ikE}\prod_{j=1}^N\left(e^{-ik\hbar\omega/2}\sum_{n_j=0}^{\infty}e^{-ik\hbar\omega n_j}\right)=\int\frac{dk}{2\pi}e^{ikE}\left(\frac{e^{-ik\hbar\omega/2}}{1-e^{-ik\hbar\omega}}\right)^N,$$
where this last step involves summing a divergent geometric series, declaring
$$\sum_{\ell=0}^{\infty}e^{-i\alpha\ell} = \frac{1}{1-e^{-i\alpha}}$$
and ignoring the fact that this series doesn't converge in a conventional sense.
This simplifies to
$$\Omega(E) = \int\frac{dk}{2\pi}e^{N(ik(E/N) - \log(2i\sin(k\hbar\omega/2)))}$$
which is solved using the saddle-point approximation. The maximum of the argument of the exponential occurs at a value
$$k_0 = \frac{1}{\hbar\omega i}\log\frac{\frac{E}{N}+\frac{\hbar\omega}{2}}{\frac{E}{N}-\frac{\hbar\omega}{2}}$$
Which is clearly imaginary, despite the fact that in a Fourier Transform $k$ is supposed to be a real number!
In spite of all this, if you evaluate the integral using the saddle point approximation at $k=k_0$, you get the same form for $\Omega(E)$ that one derives through the garden-variety combinatorial argument!
My question(s) are

Why does this work? Specifically:

*

*Why does it make sense to write the convergence of a divergent geometric series in the form given here (is this relying on some sense of convergence other than the typical one, and if so, what? And what does that imply about convergence in stat mech?), and


*Why can you use the saddle point approximation when the maximum value does not occur in the space over which you are integrating?

Answers to this question might rely on appeals to other situations in which this math occurs and has been rationalized, physically if not mathematically.
 A: I) If we expect $\Omega(E)$ to depend analytically on the variable $\hbar\omega>0$ extended to (parts of) the complex plane, then we may regularize by introducing an $i\epsilon$ prescription, and substitute 
$$\tag{1} \hbar\omega ~\longrightarrow ~ \hbar\omega (1-i\epsilon). $$
The variable 
$$\tag{2} q~:=~ e^{-i\hbar\omega k}~\longrightarrow ~ e^{-(i+\epsilon)\hbar\omega k} $$
in the geometric series
$$\tag{3}  \sqrt{q}\sum_{n=0}^{\infty} q^n $$
will then have
$$\tag{4} |q|~<~1,$$
so that the geometric series (3) is convergent. Then all steps in Schwabl's derivation of $\Omega(E)$ are mathematically well-defined. At the end of the calculation of $\Omega(E)$, we may put $\epsilon=0$.
II) Concerning a complex (as opposed to real) stationary solution in the method of steepest descent/ stationary phase method/saddle-point method, this is just part of the method. For a rigorous argument, one would have to consult the proof of the method. Heuristically, it is because when one evaluates the Gaussian integral over 'quantum fluctuations'
$$\tag{5} \int_{\mathbb{R}} \! dx~ e^{-\frac{a}{2}x^2+bx} ~=~\sqrt{\frac{2\pi}{a}}\exp\left(\frac{b^2}{2a}\right),$$
for two complex constants $a,b\in\mathbb{C}$, one only needs the condition 
$$\tag{6} {\rm Re}(a)~>~0$$
to ensure convergence of the integral (5). There is no need to also assume that the stationary solution $\frac{b}{a}$ is real. Equation (5) follows from the fact that 
$$\tag{7} \alpha\int_{\mathbb{R}} \! dx~e^{-\frac{1}{2}(\alpha x+\beta)^2}
~=~\int_{\gamma} \! d(\alpha x+\beta)~ e^{-\frac{1}{2}(\alpha x+\beta)^2}
~=~ \int_{\mathbb{R}} \! dx~ e^{-\frac{1}{2}x^2}~=~\sqrt{2\pi}$$ 
for any straight line in the complex plane 
$$\tag{8} \gamma(x)~=~\alpha x+\beta, \qquad \alpha,\beta~\in~\mathbb{C},\qquad 
x~\in~\mathbb{R}, $$ 
with slope 
$$\tag{9} |\arg(\alpha)|<\frac{\pi}{4},$$ 
because in that case, it is possible to close the contour along exponentially suppressed arcs.
