  We study the thermal partition function of level $k$ U(N) Chern-Simons
theories on $S^2$ interacting with matter in the fundamental representation. We
work in the 't Hooft limit, $N,k\to\infty$, with $\lambda = N/k$ and $\frac{T^2
V_{2}}{N}$ held fixed where $T$ is the temperature and $V_{2}$ the volume of
the sphere. An effective action proposed in arXiv:1211.4843 relates the
partition function to the expectation value of a `potential' function of the
$S^1$ holonomy in pure Chern-Simons theory; in several examples we compute the
holonomy potential as a function of $\lambda$. We use level rank duality of
pure Chern-Simons theory to demonstrate the equality of thermal partition
functions of previously conjectured dual pairs of theories as a function of the
temperature. We reduce the partition function to a matrix integral over
holonomies. The summation over flux sectors quantizes the eigenvalues of this
matrix in units of ${2\pi \over k}$ and the eigenvalue density of the holonomy
matrix is bounded from above by $\frac{1}{2 \pi \lambda}$. The corresponding
matrix integrals generically undergo two phase transitions as a function of
temperature. For several Chern-Simons matter theories we are able to exactly
solve the relevant matrix models in the low temperature phase, and determine
the phase transition temperature as a function of $\lambda$. At low
temperatures our partition function smoothly matches onto the $N$ and $\lambda$
independent free energy of a gas of non renormalized multi trace operators. We
also find an exact solution to a simple toy matrix model; the large $N$
Gross-Witten-Wadia matrix integral subject to an upper bound on eigenvalue
density.
