  We rigorously compute the integrable system for the limiting
$(N\rightarrow\infty)$ distribution function of the extreme momentum of $N$
noninteracting fermions when confined to an anharmonic trap $V(q)=q^{2n}$ for
$n\in\mathbb{Z}_{\geq 1}$ at positive temperature. More precisely, the edge
momentum statistics in the harmonic trap $n=1$ are known to obey the weak
asymmetric KPZ crossover law which is realized via the finite temperature Airy
kernel determinant or equivalently via a Painlev\'e-II integro-differential
transcendent, cf. \cite{LW,ACQ}. For general $n\geq 2$, a novel higher order
finite temperature Airy kernel has recently emerged in physics literature
\cite{DMS} and we show that the corresponding edge law in momentum space is now
governed by a distinguished Painlev\'e-II integro-differential hierarchy. Our
analysis is based on operator-valued Riemann-Hilbert techniques which produce a
Lax pair for an operator-valued Painlev\'e-II ODE system that naturally encodes
the aforementioned hierarchy. As byproduct, we establish a connection of the
integro-differential Painlev\'e-II hierarchy to a novel integro-differential
mKdV hierarchy.
