statistical mechanics - Canonical averages in a Fermi gas aka generalized Fermi-Dirac distribution

I am in the process of applying Beenakker's tunneling master equation theory of quantum dots (with some generalizations) to some problems of non-adiabatic charge pumping. As a part of this work I encounter thermal averages of single-particle quantities with a fixed total number of electrons. They are pretty straightforward to derive, as recently discussed on Physics.SE, see Combinatorial sum in a problem with a Fermi gas.

I have trouble finding references to prior work for this basic quantum statistics problem. Can you suggest some relevant references in this context?

Progress report:

Writing this stuff up, I realized that I'm basically asking for finite electron number and finite level spacing generalization of the Fermi function:

$\langle \nu_k \rangle = Z_n^{-1} \sum_{n=0}^{n-1} (-1)^{n-m} Z_m e^{-\beta \epsilon_k( n-m)}$ where $Z_n$ is the $n$-electron canonical partition function.

For either $n \gg 1$ or and $\epsilon_{k+1}-\epsilon_k \ll \beta^{-1}$ this reduces to the standard Fermi-Dirac distribution:

$\langle \nu_k \rangle= \frac{1}{1+z_0 e^{-\beta \epsilon_k}}$

The advantage of the new formula is that it only depends on $k$ via $e^{-\beta \epsilon_k}$ with all the combinatorics hidden in $Z_n$. This is unlike the textbook expression summarized in Wikipedia, or Beenakker's implicit function $F(E_p |n)$ (which is $\langle \nu_p \rangle$ in notation).

The question still stands: I don't believe the first formula above is new, what is the right reference?

Answer

A particle-number projection method which is based on Fourier transform of the grand-canonical partition function has been developed by Ormand et al. Phys. Rev. C 49, 1422 (1994) in the context of nuclear quantum Monte-Carlo simulations. This gives a closed form formula which scales quadratically in the number of levels. It has been used by Alhassid et al. Phys. Rev. B 58, R7524 (1998) to describe (numerically) deviations from Fermi-Dirac distribution, see Eq. (140) in Alahssid's review paper.

Since the approach described in the question appears to be new, I've posted it to the arXiv. For equidistant (1d harmonic oscillator) spectrum, generalization of Fermi-Dirac distribution leads to an intriguing series of polynomials and Rogers-Ramanujan partial theta function, see math.SE discussion.