The generating functions for the supersymmetric indices of the gauge theory such as superconformal index are often represented in terms of the unitary matrix integrals with double trace potential. In the limit of weak interactions between the eigenvalues, they can be approximated by the matrix models with the single-trace potential, i.e. the generalized Gross-Witten-Wadia model. In this work, the perturbative and non-perturbative aspects of the generic multi-critical unitary matrix models are studied by adopting the integrable operator formalism, and the multi-critical generalization of the Tracy-Widom distribution in the context of random partitions. We obtain the universal results for the multi-critical model in the weak and strong coupling phases. The free energy of the instanton sector in the weak coupling regime, and the genus expansion of the free energy in the strong coupling regime are explicitly computed and the universal multi-critical phase structure of the model is explored. Finally, we apply our results in concrete examples of supersymmetric indices of gauge theories.