  We consider classical shadows of fermion wavefunctions with $\eta$ particles
occupying $n$ modes. We prove that all $k$-Reduced Density Matrices (RDMs) may
be simultaneously estimated to an average variance of $\epsilon^{2}$ using at
most
$\binom{\eta}{k}\big(1-\frac{\eta-k}{n}\big)^{k}\frac{1+n}{1+n-k}/\epsilon^{2}$
measurements in random single-particle bases that conserve particle number, and
provide an estimator for any $k$-RDM with $\mathcal{O}(k^2\eta)$ classical
complexity. Our sample complexity is a super-exponential improvement over the
$\mathcal{O}(\binom{n}{k}\frac{\sqrt{k}}{\epsilon^{2}})$ scaling of prior
approaches as $n$ can be arbitrarily larger than $\eta$, which is common in
natural problems. Our method, in the worst-case of half-filling, still provides
a factor of $4^{k}$ advantage in sample complexity, and also estimates all
$\eta$-reduced density matrices, applicable to estimating overlaps with all
single Slater determinants, with at most $\mathcal{O}(\frac{1}{\epsilon^{2}})$
samples, which is additionally independent of $\eta$.
