Degenerate ground state of Hamiltonian from analytical perspective Suppose I have a Hamiltonian that depends on the continuous vector parameter $\boldsymbol{\theta}$, and the ground state corresponds to line/plane or some other $1$ to $p-1$ dimensional subspace of the $p\,$ dimensional space of my parameter vector.
And suppose I wish to find the ground state energy using a statistical mechanics/analysis method:


*

*I write my partition function $Z = \displaystyle\int d\boldsymbol{\theta}\,e^{-\beta H(\boldsymbol{\theta})}$

*I use the formula for the internal energy $U = -\dfrac{\partial}{\partial \beta} \ln(Z)$

*And I take the limit of $\beta \rightarrow \infty$ corresponding to $T\rightarrow 0$
This should give me the ground state energy. But, I am suspicious of whether the underlying mathematics suggests this is a fishy thing to do in the case of a degenerate ground state.
From my limited understanding of analysis:
I take:
$$U = \lim_{\beta \rightarrow \infty} \frac{\displaystyle \int d\boldsymbol{\theta}\,H(\theta)e^{-\beta H(\boldsymbol{\theta})}}{Z}$$
and this is handled with the saddle point method: Watson's lemma etc.
But while I understand how the saddle point method works (in terms of using a Taylor expansion of the Hamiltonian) for a single saddle point (this stack exchange post has a good overview) I don't think it works when there is a dimension $\dim \mathfrak g>0$ subspace for the ground state.
Would it be appropriate to do a Landau-type expansion of the Hamiltonian (though I have only come across Landau expansions of the free energy; and with macroscopic order parameters rather parameters on which the Hamiltonian depends!). But still, if I have a continuous region of the parameter space are ground states, then how to consider this analytical approach?
 A: *

*If the Hamiltonian is (continuously) degenerated, there is going to be a flat direction valley in the profile $\theta\mapsto H(\theta)$. Close to the bottom of the valley (=the ground state submanifold $\mathfrak{g}$), the equi-energy-surfaces/lines are locally almost aligned.

*Let us for simplicity assume that we can globally choose adapted parameter coordinates $$\theta = \{x,y\}\in X\times Y$$ such that 
$$ \frac{\partial H}{\partial y}~\equiv~0; $$
such that the ground state submanifold 
$$\mathfrak{g}~=~\{x\in X \mid  \frac{\partial H}{\partial x}=0\}\times Y $$ is parametrized by the $y$-coordinates; and such that the Hessian 
$$\left. \frac{\partial^2 H}{\partial x^2}\right|_{\mathfrak{g}}$$ is positive definite.

*Change of coordinates produce a Jacobian factor 
$$d\theta~=~dy~ dx~J(x,y),$$ 
so that OP's partition function becomes 
$$Z~:=~\int \!d\theta~e^{-\beta H(x)}~=~\int_{Y}\!dy~ \int_{X}\!dx~J(x,y)~e^{-\beta H(x)}, \qquad \beta\text{ large}.  $$ 

*The integration over the $x$-parameters corresponds to a non-degenerate problem (with some external parameters $y$), which OP already knows how to calculate using the method of steepest descent.

*There remains an ordinary integration over the $y$-parameters.  [The $y$-integration may or may not be doable in analytic form. If the Jacobian does not depend on $y$, then the $y$-integration just gives a volume factor ${\rm Vol}(Y)={\rm Vol}(\mathfrak{g})$.]
