  In the current era of precision quantum many-body physics, one of the most
scrutinized systems is the unitary limit of the nonrelativistic spin-$1/2$
Fermi gas, due to its simplicity and relevance for atomic, condensed matter,
and nuclear physics. The thermodynamics of this strongly correlated system is
determined by universal functions which, at high temperature, are governed by
universal virial coefficients $b_n$ that capture the effects of the $n$-body
system on the many-body dynamics. Currently, $b_2$ and $b_3$ are well
understood, but the situation is less clear for $b_4$, and no predictions have
been made for $b_5$. To answer these open questions, we implement a
nonperturbative analytic approach based on the Trotter-Suzuki factorization of
the imaginary-time evolution operator, using progressively finer temporal
lattice spacings. Implementing these factorizations and automated algebra
codes, we obtain the interaction-induced change $\Delta b_n$ from weak coupling
to unitarity. At unitarity, we find: $\Delta b_3 = -0.356(4)$, in agreement
with previous results; $\Delta b_4 = 0.062(2)$, in agreement with all previous
theoretical estimates but at odds with experimental determinations; and $\Delta
b_5 = 0.078(6)$, which is a prediction. We show the impact of those answers on
the density equation of state and Tan contact, and track their origin back to
their polarized and unpolarized components.
