# Double Sommerfeld Expansion

Consider the Fermi-Dirac expansion for an arbitrary function $$f(\epsilon)$$:

$$I(f)=\int_0^\infty d\epsilon\frac{f(\epsilon)}{e^{\beta(\epsilon-\mu)}+1}$$

The large $$\beta$$ expansion of this quantity is called the Sommerfeld expansion, and is well known to have the following all orders in $$1/\beta$$ expansion:

$$I(f)=\int_0^\mu f(\epsilon)d\epsilon+\sum_{n=1}^\infty\frac{1}{\beta^{2n}}(2-2^{2-2n})\zeta(2n)f^{(2n-1)}(\mu)+O(e^{-\beta \mu})$$

Let us now define the analogous 2-dimensional function

$$J(f)=\int_0^\infty d\epsilon_1\int_0^\infty d\epsilon_2 \frac{f(\epsilon_1,\epsilon_2)}{e^{\beta(\epsilon_1-\mu)}+1}\frac{1}{e^{\beta(\epsilon_2-\mu)}+1}$$

Is there a similar asymptotic expansion for large $$\beta$$ for this when $$f(\epsilon_1,\epsilon_2)\neq g(\epsilon_1)h(\epsilon_2)$$ for some functions $$g,h$$?

Of course in the factorizable case, one can just do the Sommerfeld expansion to each function $$g,h$$, seperately. But when the function does not factorize, how is this expansion done? For concreteness, we can consider the simple case

$$J(|{\epsilon_1-\epsilon_2}|)=\frac{\mu^3}{3}+\frac{\pi^2\mu}{3\beta^2}-\frac{4\zeta(3)}{\beta^3}+O(e^{-\beta\mu})$$

where I computed this by doing the integral exactly and then expanding in large $$\beta$$. How can one recover these perturbative terms from a Sommerfeld-like expansion?