In this article, we study the large $n$ asymptotic expansions of $n\times n$ Toeplitz determinants whose symbols are indicator functions of unions of arc-intervals of the unit circle. In particular, we use an Hermitian matrix model reformulation of the problem to provide a rigorous derivation of the general form of the large $n$ expansion when the symbol is an indicator function of either a single arc-interval or several arc-intervals with a discrete rotational symmetry. Moreover, we prove that the coefficients in the expansions can be reconstructed, up to some constants, from the Eynard-Orantin topological recursion applied to some explicit spectral curves. In addition, when the symbol is an indicator function of a single arc-interval, we provide the corresponding normalizing constants using a Selberg integral and illustrate the theoretical results with numeric simulations up to order $o\left(\frac{1}{n^4}\right)$. We also briefly discuss the situation when the number of arc-intervals increases with $n$, as well as more general Toeplitz determinants to which we may apply the present strategy.